Showing posts with label adaptor. Show all posts
Showing posts with label adaptor. Show all posts

Tuesday, December 03, 2013

Illumina HiSeq2000 adaptor and sequencing


I've referred the following material when making the figure:

Just one learning note:

If insert is not long enough (i.e. shorter than the read length), R1 will have contamination from Rd2 SP (e.g. AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC) and R2 will have contamination from reverse complementary of Rd1 SP (e.g. AGATCGGAAGAGCG)

So, basically you just need to check if the shared complementary part of Rd1 and Rd2 SP (which is AGATCGGAAGAGC) occurs in the reads. If yes, simply trim it and its following part (if any).

Note: if you don't understand the "shared complementary part", please refer my previous blog on Illumina adaptor. Here is the link: http://onetipperday.blogspot.com/2013/06/illumina-hiseq2000-adaptor.html

Here is one solution of howto remove adaptor contamination:
1. save the complementary part into a fasta file, e.g. adaptor.fa
>adaptor_complementary_part
AGATCGGAAGAGC
2. run fastq-mcf to remove adaptor
fastq-mcf -o filted -x 10 -l 15 -w 4 -u adaptor.fa input.fq.gz

Tuesday, June 04, 2013

Illumina HiSeq2000 adaptor

It's made just to help myself understanding the adaptor structure of Illumina sequencing. I will add the part of sequencing (on flow cell) later.

Continue to read the following part:

Monday, May 27, 2013

Unstable output from fastq-mcf (adaptor removal)

Below is what I wrote to the author regarding the problem I just found. Erik suggested to use "-t 0" to fix the problem (see info at the end)

================

I am really frustrated with a problem I am facing while using the fastq-mcf for adaptor removal -- the same read in different files with different number of reads will get different trimming result. WHAT'S THE PROBLEM?

Here is the test example:

adaptor:

$ cat adaptor.fa 
>adapter1
CAAGCAGAAGACGGCATACGAGATTCAAGTGTGACTGGAGTTCCTTGGCACCCGAGAATTCCA
>RC_adapter1
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACACTTGAATCTCGTATGCCGTCTTCTGCTTG

Input file:
$ head test.fa
@HWI-ST560:107:D25EJACXX:3:1301:11693:29455 1:N:0:GCCTAA
CAGTGCAATGTTAAAAGGGCATTTGGAATTCTCGGGTGCCAAGGAACTCC
+
CCCDFFFFHHHHGJJJJJJJJJJJJJJJJJJJJJJIHHJJIJJJIJJIIJ
@HWI-ST560:107:D25EJACXX:3:1101:1420:1954 1:N:0:TTAGGC
ATACAGTCCGACGATCTGGAATTCTCGGGTGCCAAGGAACTCCAGTCACT
+
@@@DDDDDHHG?FDGIIBEHAFGIIHFHIHGIIEHGGGHFHIIGGGGGHI

When I include different number of reads, the result are differnt for the 1st read:

$ head -n10000 test.fa > test1.fq; ~/bin/ea-utils.1.1.2-537/fastq-mcf -o test1.out -l 16 -q 15 -w 4 -x 20 -u -S -P 33 adaptor.fa test1.fq; head test1.out

@HWI-ST560:107:D25EJACXX:3:1301:11693:29455 1:N:0:GCCTAA
CAGTGCAATGTTAAAAGGGCATTTGGAATTCTCGGGTGCCAAGGAACTCC
+
CCCDFFFFHHHHGJJJJJJJJJJJJJJJJJJJJJJIHHJJIJJJIJJIIJ

This is not correct, as the yellow highlight part is from RC_adapter1 obviously. 

But if I set a smaller number of sample, it will be correct. See below for "head -n100":

$ head -n100 test.fa > test1.fq; ~/bin/ea-utils.1.1.2-537/fastq-mcf -o test1.out -l 16 -q 15 -w 4 -x 20 -u -S -P 33 adaptor.fa test1.fq; head test1.out

@HWI-ST560:107:D25EJACXX:3:1301:11693:29455 1:N:0:GCCTAA
CAGTGCAATGTTAAAAGGGCATT
+
CCCDFFFFHHHHGJJJJJJJJJJ

And if I set "-C 0", it will get correct result as well (even for a larger number of sample):

$ head -n10000 test.fa > test1.fq; ~/bin/ea-utils.1.1.2-537/fastq-mcf -o test1.out -l 16 -q 15 -w 4 -x 20 -u -S -P 33 -C 0 adaptor.fa test1.fq; head test1.out

@HWI-ST560:107:D25EJACXX:3:1301:11693:29455 1:N:0:GCCTAA
CAGTGCAATGTTAAAAGGGCATT
+
CCCDFFFFHHHHGJJJJJJJJJJ

This is very confusing! Please help me on it.

=========== reply from Erik========
fastq-mcf will not bother to clip if the adaptor is rare (< .25% presence in sampling).   The reason is that the dangers of clipping (removing valid sequence) can outweigh the benefits (removing invalid sequence) as the number of adaptor sequences get fewer and fewer.

To override this, you set -t 0 (always clip everything)

Wednesday, March 14, 2012

adaptor remover/clipper with mismatch allowed


$ head -n12 RNAseq/mouse_adult_wt_smallRNAseq_76nt_strand_ox.fastq 
@HWUSI-EAS1533_0026_FC:6:1:1204:14021#0/1
TTCCTATCAATAAATCTCATGTGANNGCAGCTCGTANNCCNTCTNCNNNTNNNNANNNNANANCNNNNNANNNNNN
+HWUSI-EAS1533_0026_FC:6:1:1204:14021#0/1
gggaegggggfdggfcaedbBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
@HWUSI-EAS1533_0026_FC:6:1:1204:19652#0/1
TGGTATAAACTCTTCTAAAAGACTNNAGCTCGTATGNNGCNTTCNGNNNGNNNNANNNNANANANNNNNANNNNNN
+HWUSI-EAS1533_0026_FC:6:1:1204:19652#0/1
c^^[^a`ddYaddadf_YafSOTZDDYZK\M`W_`BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
@HWUSI-EAS1533_0026_FC:6:1:1205:19274#0/1
TCAGAAATGTCCTGTGAAACTCGCNNAATTTCGTATGCCGCCTTCTGNNTGNAAANNNAAAACAANNNNANNNNNN
+HWUSI-EAS1533_0026_FC:6:1:1205:19274#0/1
fV^V]b]`ZSTXdRbc]dac`BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB


The basic question is to trim/clip the part after the adaptor sequence (e.g. TCGTATGCCGTCTTCTGCTTG in above example) in the FASTQ reads. It's easy if there is no mismatch. For case with mismatch, here is few solutions I got from google:
  1. FAR (http://sourceforge.net/apps/mediawiki/theflexibleadap/index.php?title=Main_Page), which uses Needleman-Wunsch alignment (global alignment) to align the adaptor sequence(s) to the read sequences and find the first best alignment.

    far -s stdin.fastq -t stdout.fq -f fastq-sanger -as TCGTATGCCGTCTTCTGCTTG --cut-off 5 --min-overlap 15 --phred-pre-trim 20 --min-readlength 20 --trim-end right --adaptive-overlap yes --max-uncalled 30 --nr-threads 8 --log-level ALL
  2.  Scythe only checks for 3'-end contaminants, up to the adapter's length into the 3'-end. For reads with contamination in any position, the program TagDust (http://genome.gsc.riken.jp/osc/english/dataresource/) is recommended. Scythe has the advantages of allowing fuzzier matching and being base quality-aware, while TagDust has the advantages of very fast matching (but allowing few mismatches, and not considering quality) and FDR. TagDust also removes contaminated reads entirely, while Scythe trims off contaminants.

  3. A possible pipeline would run FASTQ reads through Scythe, then TagDust, then a quality-based trimmer, and finally through a read quality statistics program such as qrqc (http://bioconductor.org/packages/devel/bioc/html/qrqc.html) or FASTqc (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/).
  4. trimLRPatterns() function in BioConductor. Here is an example.
  5. Another very useful toolkit is fastx_toolkit, which includes a set of tools for handling fasta and fastq format. fastx_clipper  is for clipping adaptor, but from what I checked it does not allow mismatch (e.g. N is counted as 'match' in fastx_clipper, but in FAR it's mismatch).  Running the following code on above reads only clips the last read. The first one should also be clipped. Don't know why. (Just get reply from Gordon that fastx_clipper does allow mismatch and also count 'N' as mismatch. Just the first two sequences have high error rate (=number_of_mismatch/adaptor_alignment_length, e.g. 33% in this case), so it's not counted as an adaptor. 
$ head -n12 RNAseq/mouse_adult_wt_smallRNAseq_76nt_strand_ox.fastq | fastx_clipper -a TCGTATGCCGTCTTCTGCTTG -c -n -i - 
@HWUSI-EAS1533_0026_FC:6:1:1205:19274#0/1
TCAGAAATGTCCTGTGAAACTCGCNNAATT
+HWUSI-EAS1533_0026_FC:6:1:1205:19274#0/1
fV^V]b]`ZSTXdRbc]dac`BBBBBBBBB