how to remove first 15 bases of a sequence using unix or qiime

[colinn]

So I'm pretty sure all of my sequences have a 15 base pair adaptor sequence at the beginning, and I'd like to remove it, so then my downstream analysis will play nicely within QIIME.

I've tried changing my mapping file by adding 15bp to each barcode of (NNNNNNNNNNNNNNN), but QIIME says thats not a cool thing to do.

I'm interested in other ways to do this, and again for the .qual file.

[luc]

Hi,

I am not entirely sure what you want to achieve. But is sounds like BBduk or BB reformat could easily do it (http://seqanswers.com/forums/showthread.php?t=58221 ) as could do some simple scripts.

[GenoMax]

Can you post an example 2-3 sequences from the fasta and qual files? BBduk from BBMap can manage the fasta file trim (forcetrimleft=15) for sure but I am not certain about the qual file. You can try it yourself.

We can always do an awk/perl solution if this does not work.

[Brian Bushnell]

You can do that with BBDuk:

bbduk.sh in=reads.fq out=trimmed.fq ftl=15

If you have fasta + qual files:

bbduk.sh in=reads.fa qfin=reads.qual out=trimmed.fa qfout=trimmed.qual ftl=15

If not all of the reads contain that, and you only want to trim the reads that have a specific sequence, then use the flags "ktrim=l literal=X k=15" instead, where X is the sequence you want to trim. If you have lots of sequences, you can separate them with commas.

Edit - double-ninja'd!

@GenoMax. Sorry the numbering is a little confusing for force-trim; these are the definitions:

	/** Trim left bases of the read to this position (exclusive, 0-based) */
	private final int forceTrimLeft;
	/** Trim right bases of the read after this position (exclusive, 0-based) */
	private final int forceTrimRight;
	/** Trim this many rightmost bases of the read */
	private final int forceTrimRight2;

So to make a 100bp read 85 bp long by trimming the right end, you would say "ftr=84" or "ftr2=15"; to do so by trimming the left side, it would be "ftl=15".

[GenoMax]

Brian: Shouldn't that be ftl=14 (0-based correct)?

[colinn]

Hi Geno-

Here are some sample lines of the .fna file, its pretty standard.

>SRR1502445.1
GACTACACGTAGTATACGAGTGCGTTCCTGCGCTTATTGATATGCTTAAGTTCAGCGGGTAGTCTCACCCGATTTGAGGTCAAGGTTTGTGGGTCGAGTCACAACTCGAACATCGTCTTTGTACAAAGACGGTTGGAAGCGGGTTCCAAGGCAACACAGGGGGATAGGNNNNNNNNNNNNNNNNNNNNNNN
>SRR1502445.2
GACTACACGTAGTATACGAGTGCGTTCCTGCGCTTATTGATATGCTTAAGTTCAGCGGGTAGTCTCACCCGATTTGAGGTCAAGGTTTGTGGGTCGAGTCACAACTCGAACATCGTCTTTGTACAAGACGGTTGGAAGCGGGTTCCAAGGCACACAGGGGATAGGNNN
>SRR1502445.3
GACTACACGTAGTATACGAGTGCGTTCCTGCGCTTATTGATATGCTTAAGTTCAGCGGGTAGTCTCACCCGATTTGAGGTCAAGGTTTGTGGGTCGAGTCACAACTCGAACATCGTCTTTGTACAAAGACGGTTGGAAGCGGGTTCCAAGGCACACAGGGGATAGGNNN
>SRR1502445.4
GACTACACGTAGTATACGAGTGCGTTCCTGCGCTTATTGATATGCTTAAGTTCAGCGGGTAGTCTCACCCGATTTGAGGTCAAGGTTTGTGGGTCGAGTCACAACTCGAACATCGTCTTTGTACAAAGACGGTTGGAAGCGGGTTCCAAGGCACACAGGGGATAGGNNNNNNNNNNN

Also- thanks for the BBDuk suggestions, however it'd be great to do this with just the base terminal commands rather than installing a new bioinformatics package.

[GenoMax]

No installation is required for BBMap. Download, uncompress and start using (I assume you have Java installed).

[SES]

No installation is required for BBMap. Download, uncompress and start using (I assume you have Java installed).

Those are conflicting statements, you can't accomplish not installing anything by installing something . I think I know what you mean though, no external dependencies. We don't really need anything complicated for this task, perl or awk should do fine as mentioned above.

Here is an example using Perl:

perl -0076 -e 'while (<>) { chomp; my ($name, @parts) = split /\n/; my $seq = join "", @parts; next unless defined $name && defined $seq; substr($seq, 0, 15, ""); print join "\n", ">".$name, "$seq\n" }' seqs.fas

You would repeat this for seqs.qual.

[dariober]

So I'm pretty sure all of my sequences have a 15 base pair adaptor sequence at the beginning

Rather than clipping the first 15bp no matter what is in there I would explicitly look for the adapter sequence and trim that. There are tools for that, I use cutadapt, but I'm sure the BBtools can do it as well.

[GenoMax]

Rather than clipping the first 15bp no matter what is in there I would explicitly look for the adapter sequence and trim that. There are tools for that, I use cutadapt, but I'm sure the BBtools can do it as well.

@dario: See this thread for background.

[pmiguel]

Here is an example using Perl:

perl -0076 -e 'while (<>) { chomp; my ($name, @parts) = split /\n/; my $seq = join "", @parts; next unless defined $name && defined $seq; substr($seq, 0, 15, ""); print join "\n", ">".$name, "$seq\n" }' seqs.fas

You would repeat this for seqs.qual.

Are we golfing?

perl -i'.bak' -pe 's/(>.*?\n).{15}/$1/' seqs.fas

perl -i'.bak' -pe 's/(>.*?\n)(\s*\d\s+){15}/$1/' seqs.qual

Warning: I did not test these.

--
Phillip