This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Remove site with only missing data

Hi! I have a multi-individual fasta file, and from this file I'm trying to remove all sites (i.e. nucleotides at the same position in the different fasta files) that contain only Ns. Do you know how I can do this? Thanks

fasta
• 1,439 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Leane • Leane

For example, if I have 4 individuals in my files:

>seq1
AAATTTGNNNNNNCN
>seq2
AAAATTCNNNNNGCA
>seq3
AAATTAGNNNANGCA
>seq4
AAATTACNNNNNNCA

I want to remove positions with 'N' from all my individuals' sequences, and obtain something like this:

>seq1
AAATTTGNNCN
>seq2
AAAATTCNGCA
>seq3
AAATTAGAGCA
>seq4
AAATTACNNCA

And keep my alignment intact.

• 0 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Leane • Leane

1 answer

assuming two lines per fasta record.

linearize the fasta, print one base per sequence, get the places where there is a N, count the places where there is 4 'N'. create a pattern for cut

CUT=`cat input.fa | paste - - | cut -f2 | while read S; do echo "${S}" | fold -w1 | awk '($1=="N") {print NR}' ; done  | sort | uniq -c | awk '($1==4) {print $2;}' | paste -sd','`

cat input.fa | while read N; do echo $N; read S ; echo "${S}" | cut --complement -c "${CUT}" ; done

>seq1
AAATTTGNNCN
>seq2
AAAATTCNGCA
>seq3
AAATTAGAGCA
>seq4
AAATTACNNCA
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

This script fails to detect the CUT object (which is therefore empty) and as a result, it cannot find the positions where all individuals have missing data (N). Sorry, I'm a beginner, but is it because I have long sequences that span more than one line? I've tried concatenating my sequences, and it still doesn't work. I also tried extracting the headers and sequences of the file using awk (with this command: awk 'BEGIN {RS=">"; ORS=""; FS="\n"} {if (NF>1) {print $2}}'), but with no result. Do you know how to fix my issue? Thanks in advance!

• 0 views • link
written by Leane • Leane

e I have long sequences that span more than one line?

yes, see : to convert a fasta to two-lines records

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Log in to answer this question.

More posts like this