Deleting sequences based on character frequency

Shell Programming and Scripting
1

This is what I would like to accomplish, I have an input file (file A) that consist of thousands of sequence elements with the same number of characters (length), each headed by a free text header starting with the chevron �>' character followed by the ID (all different IDs with different lenghts) and a number (prevalence). Something like this:

I need to calculate the frequency of A,G,C,T,- and N in each position. Then, I need to evaluate the first position in the first sequence and if the frequency of the character in that position does not reach 5% the entire sequence along with the ID should be removed. If the frequency is higher then I need to determine the frequency of the character in the second position so on and so forth till the entire sequence has been scanned. Then, I should do the same thing for each and every sequence in the file.
Thus, in my example above Seq ID 1 should be removed since "C" accounts for only 1% (Prev=1) in that column (the ID is not considered in the analysis). As result of this process, the output file (File B) should only contained Sequences 2, 3 and 4 since the frequency of each character in every position along the entire sequence in the three entries is higher than 5%. The output file should have the same format (FASTA) as the input file:

Any help will be greatly appreciated.

2

Hi, try this:

awk '$1==">"&&$5>=5{print; getline; print}' infile

or more cryptic:

awk '/^>/{p=($5>=5)?1:0}p' infile
3

Would it be possible to couple the pre-processing step? This is what I would like to do:

My input file (file A) contains the following information:

Now, I would like to create a file (file B) containing only the sequences with 5 or more characters but less than 18 with their corresponding ID an frequency, form my example above:

Then couple it with the code you mentioned above.
Would that be possible?
Thanks!

4

Hi, this should do it:

WHINY_USERS=1 awk '/^>/{id=$0;next}length>=5&&length<=18{F[$1]++;if (!I[$1]) I[$1]=id} 
                   END{for(i in I)printf "%s Freq %s\n%s\n",I,F,i}' infile

Combining these two:

WHINY_USERS=1 awk '/^>/{id=$0;next}length>=5&&length<=18{F[$1]++;if (!I[$1]) I[$1]=id} 
                   END{for(i in I)if (F>=5) printf "%s Freq %s\n%s\n",I,F,i}' infile
5

If we use File A.fas contain the following sequences

ANd then used your code

awk '/^>/{p=($5>=5)?1:0}p' A.fas > Unique.fas

As result I get Seq2, 3 and 4 but not Seq1. Now, if we determine the frequency of each character in Seq1, the values are always higher than 5% (34%), therefore it should also be kept. Thus, the code is not calculating the frequency of each character in each position (column) to decide whether to keep or remove the sequence. It is only considering the prevalence, if higher than 5% then the sequence will be kept.
I have also tried your last code:

$ awk '/^>/{id=$0;next}length>=5&&length<=18{F[$1]++;if (!I[$1]) I[$1]=id}END {for(i in I)printf "%s Freq %s\n%s\n",I,F,i}' A.fas > Unique.fas

But the Unique.fas file is empty.

6

I think then that I do not fully understand your requirement. I took it that the number behind prev or Freq or whatever is the frequency or prevalence of the string itself. You seem to be suggesting that the script should look for the frequency of a particular character in that string?

7

My bad, l guess I did not express myself correctly. From my previous example:

The first character in the first sequence is A, and the frequency of that sequence is 1%. However, A is also present in the same position (1) in sequence #2 and the frequency of that second sequence is 33%. Therefore, the 'global' frequency of A in the data set in position 1 is = 34%, hence it should be kept. When you do the same for the second and thrid position you can see that the global frequency for T (second character) and C (third character) from Sequence 1, are =34%, and therfore that should be enought to keep the entire sequence.
I hope I have better explain what I would like to accomplish.
Thanks one more time!

8

Hi, I think I understand what you mean, for the first question:

awk '{ for (i=1;i<=length($5);i++) F[(substr($5,i,1)),i]+=$4; A[NR]=$0; S[NR]=$5}
     END {for (j=1;j<=NR;j++){{p=1; for(i=1;i<=length(S[j]);i++) if(F[(substr(S[j],i,1)),i]<=5) p=0} if (p){print A[j]}}}' RS=">" ORS=">" infile

If your file is really big and you run out of memory then you could try this

awk 'NR==FNR { for (i=1;i<=length($5);i++) F[(substr($5,i,1)),i]+=$4; next }
     {p=1; for(i=1;i<=length($5);i++) if(F[(substr($5,i,1)),i]<=5) p=0} p' RS=">" ORS=">" infile infile

Output first Sample:

> ID 2 Prev 31
A-TGCTAGCTACGTCGTACGT
> ID 3 Prev 30
A-TGCTAGCTACGTCGTTCGT
> ID 4 Prev 30
A-TGCTAGCTACGTCGTACGA

Output second sample:

> Seq 1 Freq 1
ATC
> Seq 2 Freq 33
ACG
> Seq 3 Freq 33
TTG
> Seq 4 Freq 33
TCC
9

It is working great!

---------- Post updated at 08:33 PM ---------- Previous update was at 06:19 PM ----------

For a reason I can not understand, when I work with my actual data the freq is added to a second line instead to the firts one:

Original file

after using your code:

awk '/^>/{id=$0;next}length>=5&&length<=18{F[$1]++;if (!I[$1]) I[$1]=id}END{for(i in I)printf "%s Freq %s\n%s\n",I,F,i}' inputfile > Outputfile

I get the following:

The problem is that the FASTA format is distorted and I cannot process the data afterwards. Any ideas how to avoid that problem and instead generate the ID followed by the Freq:

Thanks!

10

Scrutinizer,

Here you have a couple of examples of the input files I am trying to process using your code

awk '{ for (i=1;i<=length($5);i++) F[(substr($5,i,1)),i]+=$4; A[NR]=$0; S[NR]=$5}
     END {for (j=1;j<=NR;j++){{p=1; for(i=1;i<=length(S[j]);i++) if(F[(substr(S[j],i,1)),i]<=5) p=0} if (p){print A[j]}}}' RS=">" ORS=">" inputfile > outputfile

Thanks once again!