2

I have a large .fna file containing millions of DNA sequences. Due to some write errors after truncating my primers, I'm left with some sequence headers with no sequence following it. An example is shown below. I need to delete these headers as well as the empty lines so I can continue with my workflow scripts.

B2S312_20116609 HWI-ST1128:217:D2HB4ACXX:8:1101:20814:45721 orig_bc=CCGCATAA new_bc=CCGCATAA bc_diffs=0
AACCATGCACCACCTGTATCTGTGTCTAACCAAAGGTCAGAACAACACAATCTCTTGTGTCCTCACAGTATGTCAAGACTT

B4S315_20116610 HWI-ST1128:217:D2HB4ACXX:8:1101:20948:45746 orig_bc=ATATGCCG new_bc=ATATGCCG bc_diffs=0
AACCATGCACCACCTGCACACGACCAACTAAATGCCACCACATCTCTGCAGTGTCGCCGTGCATGTCAAGCCTT

A3C214_20116611 HWI-ST1128:217:D2HB4ACXX:8:1101:20874:45747 orig_bc=GATCCAAC new_bc=GATCCAAC bc_diffs=0

A2C214_20116612 HWI-ST1128:217:D2HB4ACXX:8:1101:21248:45534 orig_bc=CCTAGGAT new_bc=CCTAGGAT bc_diffs=0
AACCGTGCACCACCTGTTTTCTGGCTTCCGAAGAAGAGGAACTATCTCTAGTTCTGTCCATCAATGTCAAGACCT

So in this example I would want to remove the A3C214_20116611 HWI-ST1128:217:D2HB4ACXX:8:1101:20874:45747 orig_bc=GATCCAAC new_bc=GATCCAAC bc_diffs=0 line and the empty line below it.

I've managed to only delete the empty line using sed, but cannot find an example that is suitable for my scenario.

  • sed '/bc_diffs=0\s*$/{N;d}' file.fna – Costas Jan 26 '15 at 10:39
  • To clarify, would you show the complete desired output? For both your input and output samples, you will get best results by using the format tool for code (marked with braces {}) rather than the one for quotes (marked by "). – John1024 Jan 26 '15 at 19:39
2
sed '$!N;/:.*\n$/d;P;D' <infile

The above command pulls in the Next line on every line which is !not the $last. If it encounters a pattern space in which a : colon is found and the last character is a \newline it deletes the lot. This means that for lines which contain a colon and the following line is blank both will be deleted. For all others it Prints up to the first occuring \newline in pattern space, then Deletes same before starting over with what remains. This gives you a one-line lookahead throughout the file.

  • 1
    This worked! The scripts I need to run are now running without error. Thank you. – Ting Jan 26 '15 at 18:40
0

(For the solutions below, I assume that you wanted one blank line to precede each header+sequence line. If you want it after, let me know; that is an easy change to make.)

Using sed

$ sed -rn 's/(.{106})/\n\1/p' file

B2S312_20116609 HWI-ST1128:217:D2HB4ACXX:8:1101:20814:45721 orig_bc=CCGCATAA new_bc=CCGCATAA bc_diffs=0 AACCATGCACCACCTGTATCTGTGTCTAACCAAAGGTCAGAACAACACAATCTCTTGTGTCCTCACAGTATGTCAAGACTT

B4S315_20116610 HWI-ST1128:217:D2HB4ACXX:8:1101:20948:45746 orig_bc=ATATGCCG new_bc=ATATGCCG bc_diffs=0 AACCATGCACCACCTGCACACGACCAACTAAATGCCACCACATCTCTGCAGTGTCGCCGTGCATGTCAAGCCTT

A2C214_20116612 HWI-ST1128:217:D2HB4ACXX:8:1101:21248:45534 orig_bc=CCTAGGAT new_bc=CCTAGGAT bc_diffs=0 AACCGTGCACCACCTGTTTTCTGGCTTCCGAAGAAGAGGAACTATCTCTAGTTCTGTCCATCAATGTCAAGACCT

How it works

  • -r

    -r tells sed to use extended regular expressions.

  • -n

    -n tells sed not to print any lines unless explicitly asked to.

  • s/(.{106})/\n\1/p

    If any line has at least 106 characters, we add a newline to the beginning and print it.

    Because the blank lines do not have 106 characters, they are not printed.

    Because lines without DNA sequences are shorter than 106 characters, they are not printed.

    Lines with DNA sequences (making them 106 characters or longer) are printed. Because the blank lines by themselves are not printed, we have to add a newline to the beginning of the line. Thus assures that each DNA record is preceded by a blank line.

Non-GNU sed

If using OSX (BSD), try:

sed -E -n -e 's/(.{106})/\n\1/p' file

Using awk

awk '$6 {print "\n"$0}' file

The DNA sequence is the sixth field on the line. Thus, for any line with a non-empty sixth field, we print that line preceded by a newline.

0

Just want to mention you could do this using vi:

:g/PATTERN/normal 2dd

Where PATTERN is what you're looking for. This will match PATTERN and execute 2dd (delete line - 2 times) in normal mode.

0

Using Ex / vi

Making a command of @ReulSharabin suggestion:

ex -s +'g/bc_diffs=0$/d +2' +x   file

which is similar to:

vim +':g/bc_diffs=0$/d +2' +:x   file

(+command = execute command; +x = save and quit)

Ok, ok, it is a bad idea: we should not load to memory a such huge file

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.