Ns126:Encode Methylation: Difference between revisions
Jump to navigation
Jump to search
>Shicheng |
>Shicheng |
||
Line 200: | Line 200: | ||
*haploinfo to MHB | *haploinfo to MHB | ||
../hapInfo2mld_block.pl ./Haib.merge.RD10_80up.hapinfo.txt 0.5 > Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed | ../hapInfo2mld_block.pl ./Haib.merge.RD10_80up.hapinfo.txt 0.5 > Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed | ||
Revision as of 02:57, 10 March 2016
RRBS Data Analysis to Encode Project (Fastq)
Aim 1: Haib39bioChain
- methylation haplotype block (MHB) calling with RRBS dataset
Background
- Basic:Single-end 40bp reads
- Encode: http://genome.ucsc.edu/ENCODE/downloads.html
- Encode|RRBS1: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/
- Encode|RRBS2: http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=437674359_aUhx08DjchWwtBjyCYv61EB7Yy8S&c=chr1&g=wgEncodeHaibMethylRrbs
- Encode|Methy450K: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethyl450/
- Encode|software: https://www.encodeproject.org/software/
Method
- Never Download Encode Data from UCSC. Encode Mainpage is great: https://www.encodeproject.org/search/?type=Experiment
- Download 39 biochain RRBS dataset from Encode Project (Fastq). Alignment with Bismark and merge all the BAM files
- make the haplotype calling with previous perl script. Finally MHB calling were down with Dr. Zhang perl script.
- Methyfreq based MHB calling were conducted with MethBed files download from UCSC.
Result
Summary Excel Haib Dataset manifest
- Alignment: Maybe walltime time is too short , not all the samples were aligned completely, so I extend the walltime to 72 hours.
#PBS -q glean #PBS -l nodes=1:ppn=8 #PBS -l walltime=72:00:00 bismark --bowtie2 --phred64-quals --fastq -L 30 -N 1 --multicore 2 /home/shg047/db/hg19/meth/bismark ../fastq_trim/HOT197_trimmed.fq.gz -o ../bam2
- With above setting, the alignment time usage is about 16 hours
- The time is propotional to the size of the fastq file
Sample | N(reads) | N(mapped) | P(mapping) | N(C) | N(MCPG) | N(MCHG) | N(MCHH) | N(UCPG) | N(UCHG) | N(UCHH) | P(MCPG) | P(MCHG) | P(MCHH) |
ENCFF000LVN | 42665329 | 31372051 | 73.50% | 227374577 | 9979920 | 566174 | 2568709 | 16516656 | 50531408 | 147211710 | 37.70% | 1.10% | 1.70% |
ENCFF000LWL | 37938401 | 25076517 | 66.10% | 232580318 | 11488185 | 403775 | 1382726 | 44058434 | 51038500 | 124208698 | 20.70% | 0.80% | 1.10% |
ENCFF000LVR | 36779630 | 25226017 | 68.60% | 216621540 | 12232832 | 438455 | 1846136 | 33041219 | 46480917 | 122581981 | 27.00% | 0.90% | 1.50% |
ENCFF000LVW | 36169867 | 22446797 | 62.10% | 224690704 | 12356273 | 372366 | 871608 | 41484525 | 51598506 | 118007426 | 22.90% | 0.70% | 0.70% |
ENCFF000LVB | 35375019 | 22874896 | 64.70% | 211937214 | 13517127 | 359167 | 1121727 | 31366464 | 47860817 | 117711912 | 30.10% | 0.70% | 0.90% |
ENCFF000LUP | 34214415 | 20212548 | 59.10% | 221862905 | 16119238 | 340947 | 580628 | 49637185 | 51949984 | 103234923 | 24.50% | 0.70% | 0.60% |
ENCFF000LWA | 33924878 | 20895342 | 61.60% | 226298397 | 14274674 | 319602 | 583229 | 52752860 | 53011846 | 105356186 | 21.30% | 0.60% | 0.60% |
ENCFF000LUV | 33798173 | 23173663 | 68.60% | 249630261 | 10489904 | 225941 | 428583 | 62372541 | 57063134 | 119050158 | 14.40% | 0.40% | 0.40% |
ENCFF000LWY | 33715116 | 19999502 | 59.30% | 213287381 | 13105506 | 316152 | 583822 | 46433586 | 48748868 | 104099447 | 22.00% | 0.60% | 0.60% |
ENCFF000LWW | 33642793 | 20282649 | 60.30% | 216267059 | 15823449 | 340126 | 638093 | 43292443 | 49869268 | 106303680 | 26.80% | 0.70% | 0.60% |
ENCFF000LWP | 32858930 | 20654452 | 62.90% | 219544115 | 13223974 | 319601 | 599593 | 48455363 | 50129553 | 106816031 | 21.40% | 0.60% | 0.60% |
ENCFF000LUU | 32815340 | 17866807 | 54.40% | 196534662 | 12582324 | 421909 | 799252 | 48097172 | 45287260 | 89346745 | 20.70% | 0.90% | 0.90% |
ENCFF000LVF | 31276946 | 22130843 | 70.80% | 241339808 | 11057578 | 348695 | 664774 | 64178792 | 55276052 | 109813917 | 14.70% | 0.60% | 0.60% |
ENCFF000LXB | 29681695 | 14796658 | 49.90% | 160388755 | 8182288 | 279840 | 480680 | 41788931 | 37024968 | 72632048 | 16.40% | 0.80% | 0.70% |
ENCFF000LWK | 29349446 | 12860133 | 43.80% | 137008343 | 8665157 | 302283 | 567461 | 29580229 | 31441673 | 66451540 | 22.70% | 1.00% | 0.80% |
ENCFF000LWE | 29211774 | 11162094 | 38.20% | 115410335 | 8071486 | 277005 | 511615 | 23863725 | 26859483 | 55827021 | 25.30% | 1.00% | 0.90% |
ENCFF000LVK | 27335110 | 16938233 | 62.00% | 193589572 | 12327753 | 322954 | 495345 | 51995628 | 43692466 | 84755426 | 19.20% | 0.70% | 0.60% |
ENCFF000LVO | 26190075 | 19276422 | 73.60% | 150096025 | 7012551 | 300245 | 1332640 | 19959469 | 31819734 | 89671386 | 26.00% | 0.90% | 1.50% |
ENCFF000LVA | 25660341 | 16009728 | 62.40% | 152217865 | 12748317 | 546338 | 1806399 | 26835875 | 36223859 | 74057077 | 32.20% | 1.50% | 2.40% |
ENCFF000LWD | 25467744 | 15459206 | 60.70% | 149205383 | 11335949 | 535669 | 1829440 | 30226946 | 34383379 | 70894000 | 27.30% | 1.50% | 2.50% |
ENCFF000LVU | 23511285 | 15078397 | 64.10% | 141659239 | 8631088 | 272002 | 734661 | 21915485 | 33328412 | 76777591 | 28.30% | 0.80% | 0.90% |
ENCFF000LWO | 22844877 | 14206467 | 62.20% | 153005363 | 9813445 | 247506 | 461786 | 33745829 | 35223721 | 73513076 | 22.50% | 0.70% | 0.60% |
ENCFF000LVI | 22656883 | 13335953 | 58.90% | 143634615 | 13714133 | 389619 | 650547 | 32241351 | 34451974 | 62186991 | 29.80% | 1.10% | 1.00% |
ENCFF000LUQ | 22247066 | 13650824 | 61.40% | 150112632 | 9777330 | 287617 | 504659 | 38167142 | 34898963 | 66476921 | 20.40% | 0.80% | 0.80% |
ENCFF000LUT | 22240097 | 14849726 | 66.80% | 162779228 | 9543204 | 285378 | 485629 | 41824570 | 38682843 | 71957604 | 18.60% | 0.70% | 0.70% |
ENCFF000LWH | 21943620 | 12872848 | 58.70% | 137662143 | 11585659 | 414164 | 781429 | 29526002 | 32473470 | 62881419 | 28.20% | 1.30% | 1.20% |
ENCFF000LVE | 21473578 | 14519161 | 67.60% | 157450743 | 7444292 | 221900 | 417095 | 40943018 | 36372568 | 72051870 | 15.40% | 0.60% | 0.60% |
ENCFF000LWX | 21208447 | 13218301 | 62.30% | 145215571 | 9668358 | 223882 | 356827 | 34200783 | 33840207 | 66925514 | 22.00% | 0.70% | 0.50% |
ENCFF000LVV | 21080863 | 13197949 | 62.60% | 142640926 | 10137634 | 255425 | 440789 | 32877645 | 33833167 | 65096266 | 23.60% | 0.70% | 0.70% |
ENCFF000MLE | 20338545 | 12207095 | 60.00% | 110465581 | 11201397 | 457519 | 1827464 | 19226614 | 25845153 | 51907434 | 36.80% | 1.70% | 3.40% |
ENCFF000LWS | 20241908 | 11560661 | 57.10% | 126245157 | 11387405 | 245297 | 424016 | 29669508 | 30429870 | 54089061 | 27.70% | 0.80% | 0.80% |
ENCFF000LVZ | 20058311 | 11004008 | 54.90% | 120506134 | 9791631 | 242794 | 448332 | 25470114 | 28547035 | 56006228 | 27.80% | 0.80% | 0.80% |
ENCFF000LWT | 19407909 | 10887950 | 56.10% | 119134111 | 9343415 | 203527 | 331921 | 25331949 | 28315688 | 55607611 | 26.90% | 0.70% | 0.60% |
ENCFF000MLD | 19184685 | 14418868 | 75.20% | 162604971 | 6716952 | 305536 | 587608 | 54429999 | 38467374 | 62097502 | 11.00% | 0.80% | 0.90% |
ENCFF000MLP | 18102015 | 10032177 | 55.40% | 109233809 | 10521402 | 215592 | 370895 | 25497472 | 26468867 | 46159581 | 29.20% | 0.80% | 0.80% |
ENCFF000LVJ | 17856591 | 11324104 | 63.40% | 102317536 | 8591993 | 307833 | 1097735 | 17653908 | 24511494 | 50154573 | 32.70% | 1.20% | 2.10% |
ENCFF000MLJ | 16177211 | 8675428 | 53.60% | 94678187 | 9430914 | 200777 | 318563 | 22135878 | 23265164 | 39326891 | 29.90% | 0.90% | 0.80% |
ENCFF000LUN | 10913594 | 3643074 | 33.40% | 40450477 | 3011262 | 76408 | 126647 | 9141529 | 9549119 | 18545512 | 24.80% | 0.80% | 0.70% |
ENCFF000MLM | 5127779 | 2931330 | 57.20% | 36675584 | 1281379 | 40172 | 75178 | 10725366 | 8687156 | 15866333 | 10.70% | 0.50% | 0.50% |
- MHB regions identification
#!/bin/csh #PBS -N bam2MHB #PBS -q pdafm #PBS -l nodes=1:ppn=16 #PBS -l walltime=72:00:00 #PBS -o bam2MHB.log #PBS -e bam2MHB.err #PBS -V #PBS -M shihcheng.guo@gmail.com #PBS -m abe #PBS -A k4zhang-group cd /home/shg047/oasis/Haib/sortBam # samtools cat -h header.sam -o haib.merge.bam *sort.bam samtools sort -@ 16 haib.encode.merge.bam -o haib.merge.sort.bam samtools index haib.merge.sort.bam bedtools genomecov -bg -split -ibam haib.merge.sort.bam > haib.merge.bam.pool.bed awk '$4>9 { print $1"\t"$2"\t"$3}' haib.merge.bam.pool.bed | bedtools merge -d 10 -i - > haib.RD10.genomecov.bed awk '$3-$2>80 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' haib.RD10.genomecov.bed > haib.RD10_80up.genomecov.bed
- Statistic: haib.RD10_80up.genomecov.bed
cat haib.RD10_80up.genomecov.bed|awk '{sum+=$4} END { print "N = ", NR, "Sum = ", sum, " Average = ",sum/NR}' N = 120994 Sum = 17702479 Average = 146.309
- haploinfo to MHB
cd /home/shg047/oasis/Haib/sortBam /home/shg047/oasis/Haib/mhb/haib.RD10_80up.genomecov.bed /home/shg047/oasis/Haib/mhb/haib.merge.sort.bam /home/shg047/oasis/Haib/hapInfo2mld_blocks.pl ../mergedBam2hapInfo.pl ./haib.RD10_80up.genomecov.bed haib.merge.sort.bam > Haib.merge.RD10_80up.hapinfo.txt # get hapinfo ../hapInfo2mld_block.pl ./Haib.merge.RD10_80up.hapinfo.txt 0.5 > Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed
- MHB identified with different threshold: R-square from 0.1-0.9
/home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.1.bed /home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.2.bed /home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.3.bed /home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.4.bed /home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed /home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.6.bed /home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.7.bed
R-square threshold | MHB counts |
0.1 | 14933 |
0.2 | 13367 |
0.3 | 11667 |
0.4 | 9754 |
0.5 | 8155 |
0.6 | 7683 |
0.7 | 7445 |
bedtools intersect -wa -u /home/shg047/oasis/Haib/mhb/Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed
- MHB calling based on RRBS Haib biochain data
cd /home/shg047/oasis/Haib/sortBam samtools cat -h header.sam -o haib.encode.merge.bam *sort.bam samtools sort haib.encode.merge.bam -o haib.encode.merge.sort.bam samtools index haib.encode.merge.sort.bam bedtools genomecov -bg -split -ibam haib.encode.merge.sort.bam > haib.encode.merge.bam.pool.bed awk '$4>9 { print $1"\t"$2"\t"$3}' haib.encode.merge.bam.pool.bed | bedtools merge -d 10 -i - > haib.encode.RD10.genomecov.bed awk '$3-$2>80 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' haib.encode.RD10.genomecov.bed > haib.encode.RD10_80up.genomecov.bed
- haploinfo to MHB
../hapInfo2mld_block.pl ./Haib.merge.RD10_80up.hapinfo.txt 0.5 > Haib.merge_RD10_80up.mld_blocks_r2-0.5.bed