This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Parsing large sequence file into per-sample files after barcodes/primers removed.

I have a sequence file, called seqs.fna. The sequences have had sample names assigned, barcodes and forward primers removed, and have been quality filtered. These are forward reads only, the reverse read has been discarded (don't ask). The head of the sequence file looks like this:

>SCBI.006.41.M.29.37.20140618_0 M02149:248:000000000-AHBR7:1:1101:16494:1761 1:N:0:0 orig_bc=GACTGCATCTTA new_bc=GACTGCATCTTA bc_diffs=0
GTTGGTGAACCTGCGGAAGGATCATTAAACCTTATACTTCGTGGTGCCGTCGTTGAATAGCGGTTCCTTTGTGGTTCAGCTTTTCTCCCCGCCCTGTCTGAATATTACCCATGCTTTTTGCGTACCCATTGTTTCCTCGGCAGGCTTGCCTGCCTGCTTTTCCCCTTTTTAACCTTTTGCTTTTGCAGTCAGCGTCTGAAAAAATCTTATTTGTTACATCTTTCTACAACGGTTCTCTTGGTTCTGGCTTC
>OSBS.001.23.M.29.2.20140701_1 M02149:248:000000000-AHBR7:1:1101:16261:1795 1:N:0:0 orig_bc=CTCAATGACTCA new_bc=CTCAATGACTCA bc_diffs=0
GTAGGTGAACCTGCGGAAGGATCATTAGTGAACGCCCTCACGGGCTTATACTATTCCAAACCTCTGTGTACCGTGCCCTTCGGGGCTATTTTACAAACATGGTGTAACGAACGTCATATATCATAACAAAACAAAACTTTCAACAACGGATCTCTTGGCTCTCGCATCGATGAAGAACGCAGCCGCTCAATGACTCTTTCTCGTATGCCGTCTTCTGCTTGTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
>JERC.002.M.11.0.20140805_2 M02149:248:000000000-AHBR7:1:1101:15634:1801 1:N:0:0 orig_bc=TGATCAGAAGAG new_bc=TGATCAGAAGAG bc_diffs=0
GTAGGTGAACCTGCGGAAGGATCATTAACAAATCTTGTAGTTCGCTCATCCGCGTCTTCAACAAATCCTTCCTCACCTTGTGCAATATTGACATCGTGGCAAGTTCTCTCGGGAACTTGCTGGTCTTGTGGGGGGAAAAACCTATAGGATAGGCAGCAATGCTTAACCTTTAAGCCCCCCAAAGTATTTAACCTATACGTTTTTTGTTTGTCTTAATCTTTTTATTGATACAATATATACAACTTTCATCT
>CPER.001.M.6.12.20140415_3 M02149:248:000000000-AHBR7:1:1101:16004:1804 1:N:0:0 orig_bc=ATGTCGAGAGAA new_bc=ATGTCGAGAGAA bc_diffs=0
GTAGGTGAACCTGCGGAAGGATCATTATTGAATAAACTTGGTTGGGTTGTTGCTGGCTCTTAGGAGCTTGTGCACGCTCACTCCATTTTTAACCACCTGTGCACCTTGTGTAGATCTGGACAACTCTCGGGGTAACTCGGTCAGAGAATTGCTGTGCGTAAAAGCCAGCTTTTCTTGAGTCCAGATCTATGTCTCTATATACTCTATAAGAATGTATTAGAATGTAATTGATAGGCGTTACTGCCTTTATT
>UNDE.008.21.4.7.20140707.M_4 M02149:248:000000000-AHBR7:1:1101:17724:1813 1:N:0:0 orig_bc=CCGGGATTTCAA new_bc=CCGGGATTTCAA bc_diffs=0
GTAGGTGGACCTGCGGAAGGATCATTATCGTACCACAGTGGTGCTCGGGTTGTCGCTGACCTCTTGGTCGTGCACGCCCGTGTGCTCTCAATCCATTTCACCCTTTGTGCATCACCGCGCGGGGTCTCTTCCTCTTGGCTTGCTTCAAGAGGGGGTGGTTCGCGTTTTTCATACAAACACCCTTCTAGTTTTAGTATGTCATTCTTTTGCGTTCTTACGCTATCAATACAACTTTCAACAACGGATCTCTT

The first read has the name: SCBI.006.41.M.29.37.20140618_0. The trailing _0 indicates its the first (or zero-th) sequence in the file. There are other sequences in the file with the name SCBI.006.41.M.29.37.20140618, and a trailing _N, where N is a number. Same is true for OSBS.001.23.M.29.2.20140701, and the rest of the sequence names.

I would like to separate this file into per-sample sequence files (my analysis does not necessarilly require all of these samples). Is there a straightforward way to do this?

dereplicate
• 2,593 views • link
updated by Biostar Biostar • written by caverill • caverill

I would like to separate this file into per-sample sequence files

that does not jive with the title of this post.

De-replicating sequence file after barcodes/primers removed.

So do you want to de-duplicate the sequences or separate them into files based on a per sample basis? Depending on answer you will need to change your title or text in post.

• 1 views • link
written by GenoMax GenoMax

Apologies, title is updated.

• 1 views • link
written by caverill • caverill

No worries. Just wanted to make sure what the requirement was since those two require different paths.

• 1 views • link
written by GenoMax GenoMax

What does the 'orig_bc=*' identifier stand for?

• 1 views • link
written by pfs pfs

I'm guessing original barcode, but to be honest, I do not know. Dealing with what I have, not what I want.

• 1 views • link
written by caverill • caverill

2 answers

That does look like the original illumina barcode.

It appears that this dataset is from an intermediate step in qiime pipeline. Is that correct? qiime already has a script available to do this. See here.

• 1 views • link
written by GenoMax GenoMax

That did the trick! Great catch on recognizing the qiime formatting.

• 1 views • link
written by caverill • caverill

An awk solution:

$ awk '$0 ~ "^>" {match($1, /^>([^_]+)_/, id); filename=id[1]} {print >> filename}' seqs.fna

fin swimmer

• 1 views • link
written by finswimmer finswimmer

Log in to answer this question.

More posts like this