Dinh/Dinh 2013/NOTES/2013-6-24: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
mNo edit summary
>Dinh
mNo edit summary
Line 104: Line 104:
* The long reads were not compatible with our pipeline using Bowtie2, because the aligner suppresses read name lines which are longer than 256 characters. This causes truncated original reads stored in the read name line.
* The long reads were not compatible with our pipeline using Bowtie2, because the aligner suppresses read name lines which are longer than 256 characters. This causes truncated original reads stored in the read name line.
* Added code to split long reads into multiple short reads to use Bowtie2
* Added code to split long reads into multiple short reads to use Bowtie2
* [[File:smartBisReadMapper.txt]]
* Mapping pipeline:
* Mapping pipeline:
  <nowiki>
  <nowiki>

Revision as of 01:24, 25 June 2013

Analysis of Blueprint MiSeq Test Run

Pre-processing

Standard trimming + UMI

  • Obtain UMI from the first 10bp of read 1, label both reads with UMI
  • Trim 27 bp from 5 prime end of read 1 and read 2

Adapter removal with fastq-mcf

  • Remove adapter sequences using fastq-mcf ?
  • Since the insert = (target_size= 200-280) + (2 arms ~ 64 bp) + (UMI = 10 bp) = 274-354bp, we would not have sequenced the adapters with just 250bp from each end.
  • Try removing adapters because mapping rate was ~40%.
  • Adapters:
>Linker
TTGGAGGCTCATCGTTCCTATTCAGGCAGATGTTATCGAGGTCCGAC
>Linker_rev
GTCGGACCTCGATAACATCTGCCTGAATAGGAACGATGAGCCTCCAA

fastq-mcf results

  • Note that only R1 was end-trimmed using "Linker" and only R2 was end-trimmed using "Linker_rev"
index Total reads Clipped 'end' reads % clipped Too short after clip % too short
1, R1 546,091 265,765 48.67% 16,468 3.02%
1, R2 546,091 330,822 60.58% 16,112 2.95%
2, R1 1,201,482 599,205 49.87% 46,756 3.89%
2, R2 1,201,482 733,958 61.09% 45,880 3.82%
3, R1 2,107,061 1,159,863 55.05% 30,221 1.43%
3, R2 2,107,061 1,348,459 64.00% 29,076 1.38%
4, R1 1,576,107 821,955 52.15% 16,478 1.05%
4, R2 1,576,107 1,000,442 63.48% 16,703 1.06%
5, R1 2,416,376 1,172,370 48.52% 30,785 1.27%
5, R2 2,416,376 1,470,535 60.86% 29,880 1.24%
6, R1 2,025,366 1,044,010 51.55% 17,056 0.84%
6, R2 2,025,366 1,266,104 62.51% 16,306 0.81%
7, R1 1,632,245 737,076 45.16% 12,924 0.79%
7, R2 1,632,245 398,506 24.41% 12,788 0.78%
8, R1 2,699,273 1,252,029 46.38% 24,410 0.90%
8, R2 2,699,273 668,800 24.78% 22,635 0.84%

Merge R1 and R2 with COPE

  • Create kmer_table for COPE
  • Use COPE to combine overlapping read 1 and read 2
for f in 1Index1_S1 1Index2_S2 1Index3_S3 2Index4_S6 1Index5_S4 2Index6_S7 2Index7_S8 1Index8_S5
do

	./getUMI.pl 130620_MiSeq/TES1-${f}_L001_R1_001.fastq 130620_MiSeq/TES1-${f}_L001_R2_001.fastq $f
	~/softwares/cope-src-v1.1.3/src/cope -a $f.R1.fq -b $f.R2.fq -o $f.fq -2 $f.leftR1.fq -3 $f.leftR2.fq -m 1  -t kmer_table.freq.cz -f kmer_table.freq.cz.len >cope.$f.log 2>cope.$f.error
	rm $f.R1.fq $f.R2.fq
done

COPE results

  • Without adapter removal:
index total_pairs connected_pairs connect_ratio(%) low_quality_pairs low_quality_ratio(%)
1 546,091 37,624 6.88969 239413 43.8412
2 1,201,482 76,332 6.35315 603335 50.2159
3 2,107,061 120,889 5.73733 1155859 54.8565
5 2,416,376 157,448 6.51587 1124046 46.5178
8 2,699,273 191,332 7.08828 1138412 42.1748
4 1,576,107 100,314 6.36467 805229 51.0897
6 2,025,366 131,396 6.48752 990421 48.9008
7 546,091 35,800 6.55568 239459 43.8497
  • With adapter removal:

Alignment with Bowtie2

  • The long reads were not compatible with our pipeline using Bowtie2, because the aligner suppresses read name lines which are longer than 256 characters. This causes truncated original reads stored in the read name line.
  • Added code to split long reads into multiple short reads to use Bowtie2
  • File:SmartBisReadMapper.txt
  • Mapping pipeline:
cur_dir="/media/3TB_Dinh/Blueprint"
reads_dir="/media/3TB_Dinh/Blueprint"
email="diep.hue.dinh@gmail.com"

bisReadMapper="/home/ddiep/scripts/MethylationPipeline/scripts/smartBisReadMapper.pl"
template_fwd="/media/2TB_storeA/BisRef/bisHg19/hg19.fa.bis.fwd.bowtie2"
template_rev="/media/2TB_storeA/BisRef/bisHg19/hg19.fa.bis.rev.bowtie2"
template_fa="/media/2TB_storeA/BisRef/bisHg19/hg19.fa"
soap="/home/ddiep/softwares/soap2.21release/soap"
bowtie="bowtie2"
cpg="/media/2TB_storeA/BisRef/bisHg19/C_Pos/hg19.fa.cpg.positions.txt"

INDX="1Index1_S1 1Index2_S2 1Index3_S3 2Index4_S6 1Index5_S4 2Index6_S7 2Index7_S8 1Index8_S5"

# make sure the qual_base variable is set correctly
cd $cur_dir

for s in ${INDX}
do

f="$s.fq"
g="$s.leftR1.fq"
h="$s.leftR2.fq"
n="$s-MS"
mkdir $n
echo "cd $cur_dir/$n" > $n.job
echo "$bisReadMapper -r $reads_dir/$f -W $template_fwd -C $template_rev -g $template_fa -a $bowtie -p 16 -b 33 -n $n.f1 -q 20 -l $cpg > $n.f1.statusMbias 2>$n.f1.err" >> $n.job
echo "$bisReadMapper -r $reads_dir/$g -W $template_fwd -C $template_rev -g $template_fa -a $bowtie -p 16 -b 33 -n $n.f2 -q 20 -l $cpg > $n.f2.statusMbias 2>$n.f2.err" >> $n.job
echo "$bisReadMapper -r $reads_dir/$h -W $template_fwd -C $template_rev -g $template_fa -a $bowtie -p 16 -b 33 -n $n.f3 -q 20 -l $cpg > $n.f3.statusMbias 2>$n.f3.err" >> $n.job
echo "rm *encoded" >> $n.job

done

Mapping statistics

  • The initial mapping rates (without clipping adapters were ~40%)
index read file total bases total mapped bases %mapped bases
1 merged 16078514 9618453 59.82%
1 leftover R1 104318024 50748930 48.65%
1 leftover R2 99998101 48392445 48.39%
2 merged 32606092 19767660 60.63%
2 leftover R1 230586489 105992213 45.97%
2 leftover R2 215612917 97040101 45.01%
3 merged 51734457 30624957 59.20%
3 leftover R1 405686018 160808089 39.64%
3 leftover R2 372674273 145064585 38.93%
4 merged 42833085 24714027 57.70%
4 leftover R1 300568658 131046838 43.60%
4 leftover R2 281078201 120619977 42.91%

Remove clonal reads

  • Use prof. zhang's code to remove clonal reads using UMI