Dinh/Dinh 2013/NOTES/2013-6-24: Difference between revisions
Jump to navigation
Jump to search
>Dinh mNo edit summary |
>Dinh mNo edit summary |
||
Line 41: | Line 41: | ||
| 7||546,091||35,800||6.55568||239459||43.8497 | | 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. | |||
* Added code to split long reads into multiple short reads to use Bowtie2 | |||
* Ran mapping pipeline: | |||
<nowiki> | |||
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 | |||
#sh $n.job | |||
#qsub -q home-k4zhang $n.job | |||
done | |||
</nowiki> |
Revision as of 20:56, 24 June 2013
Analysis of Blueprint MiSeq Test Run
- 2 X 250bp MiSeq PE run, total = 13,117,847 PE reads
- Noi's labnote page: http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2013-6-6
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.
- Added code to split long reads into multiple short reads to use Bowtie2
- Ran 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 #sh $n.job #qsub -q home-k4zhang $n.job done