-2

I have a fasta file (not in right format) that contains hundreds of thousands of different lengths of DNA sequences like this:

>NODE_213384_length_62_cov_8686_ID_2134025ATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAACGGAATCGAATGG>NODE_213385_length_62_cov_7933_ID_2134027ATCATCATCGAATGGAATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAAC>NODE_213386_length_62_cov_7184_ID_2134029AATGATTATTCCATTCGAGTCCATTCGATGATTCCATTCGATTCCATTCGATGATGATTGCA>NODE_213387_length_62_cov_8639_ID_2134031CAGAGCAGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAG>NODE_213388_length_62_cov_6833_ID_2134033AGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAGGTCAAT

I'd like to use a simple Linux command to extract those sequences longer than 1000bp and output in a right fasta format like this:

>NODE_213384_length_62_cov_8686_ID_2134025
ATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAACGGAATCGAATGG
>NODE_213385_length_62_cov_7933_ID_2134027
ATCATCATCGAATGGAATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAAC
>NODE_213386_length_62_cov_7184_ID_2134029
AATGATTATTCCATTCGAGTCCATTCGATGATTCCATTCGATTCCATTCGATGATGATTGCA
>NODE_213387_length_62_cov_8639_ID_2134031
CAGAGCAGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAG

Appreciate anyone who can help with this.

  • 3
    Why, in the name of all that is cretinous, did you replace the text output with an image of the same text? – jasonwryan Mar 14 '16 at 21:58
  • Yes. I replaced the original texts with images since the text did not display the ">" in the beginning of each line. – Peter Liu Mar 15 '16 at 1:48
  • wouldn't it have been simpler (and much, much better) to simply edit the question and add the missing > characters? You've now replaced usable, copy-pastable, editable, searchable text with a useless graphic. that you can't do anything with (and will be illegible on high-resolution displays) – cas Mar 15 '16 at 2:39
  • i've reverted the edit so it's text again and reformatted as code (using the {} icon in the SE editor), which preserves the > characters (which otherwise would be interpreted as block-quote) – cas Mar 15 '16 at 2:41
  • Is the title for each sequence always the same number of characters? This is unix/linux, so explain what bp are. In any case, does sequence only include ATGC, or are there ambiguous nucleotides? – Sparhawk Mar 15 '16 at 2:51
0

It looks like you have a very long line. You might run into problems with sed and awk as they eat up available memory to handle this task (depends on your file size/line length of course). So, taking the two-stage approach, but with memory constraints, you can use tr followed by one of the awk, perl, or sed solutions above.

head -20 inputfile | 
tr '>' '\n'  > stage1
perl -ne 'print ">$1 $2\n" if /^(.*?)([ACGTU]+)$/ && length($2)>1000' < stage1 > output

Once you're satisfied it works with the first 20 lines, do it for real in a single pipeline:

tr '>' '\n' inputfile | 
perl -ne 'print ">$1 $2\n" if /^(.*?)([ACGTU]+)$/ && length($2)>1000' > output

My perl script might not be as efficient as others, but should get the job done. I wrote it for clarity. Print out the label followed by a space and then the associated base pairs and newline, if and only if there is such a line that contains more than 1000 base pairs.

0

A "simple Linux command" might look something like:

sed 's/\(>[^ACGT]*_[0-9]\+\)\([ACGT]\+\)/\1\n\2\n/g' yourdnafile |egrep -B1 '^[ACGT]{1000}'

The sed part breaks down the into 2 lines per set, and the grep shows and match over 1000 and the line preceding it (-B1).

or this might be even more simple:

sed 's/\(>[^>ACGT]*\)\([ACGT]\{1000\}[ACGT]*\)/\1\n\2\n/g;s/>[^>ACGT\n]*[ACGT]\+//g' yourdnafile
  • Thanks Mel. I just tried your both commands. The first did not work. The second seemed working but no output. Since I have limited experience of Linux, could you fix it? – Peter Liu Mar 15 '16 at 14:00
  • They are both looking for 1000 characters matching A,C,G, or T. Not being familiar with DNA, perhaps I have misunderstood what 1000bp should look like? – Mel Mar 15 '16 at 14:13
0

You could do it in a two step process; first split up into the target format, then print the line pairs where the DNA sequence is long enough. E.g., on the assumption that 1000bp refers to the DNA sequence character length (rather than the value following the word "length"), it could be the following

cat inputfile | sed -e 's/>/\n>/g' -e 's/\(ID_[0-9]*\)/\n\1/g' |\
awk -F_ '$1=="NODE" { n=$0; next; } length($0)>1000 { print n; print ">" $0;}' \
> outputfile
  • Thanks Ralph. I am not a professional Linux person. Could you clarify more for where the input file should be? I guess you print the output on screen. Can you put the output in your command? – Peter Liu Mar 15 '16 at 15:11
  • Sure; I added input and output files to the command. It's broken up into three lines for readability. You can type it on a single line by not typing "backslash newline". – Ralph Rönnquist Mar 15 '16 at 20:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.