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 84: Line 84:


===Mapping statistics===
===Mapping statistics===
* The initial mapping rates (without clipping adapters were ~40%)
{| {{table}} border=1
| align="center" style="background:#f0f0f0;"|'''index'''
| align="center" style="background:#f0f0f0;"|'''read file'''
| align="center" style="background:#f0f0f0;"|'''total bases'''
| align="center" style="background:#f0f0f0;"|'''total mapped bases'''
| align="center" style="background:#f0f0f0;"|'''%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%
|-
|}
* 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.
==Remove clonal reads==
==Remove clonal reads==
* Use prof. zhang's code to remove clonal reads using UMI
* Use prof. zhang's code to remove clonal reads using UMI

Revision as of 00:58, 25 June 2013

Analysis of Blueprint MiSeq Test Run

Pre-processing

  • 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
  • Remove adapter sequences using fastq-mcf


  • 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

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

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
  • 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%
  • 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.

Remove clonal reads

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