Dinh/Dinh 2013/NOTES/2013-6-24

From ZhangLabWiki
Revision as of 22:59, 24 June 2013 by >Dinh
Jump to navigation Jump to search

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

Remove clonal reads

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