Ns126:Encode Methylation
Jump to navigation
Jump to search
RRBS Data Analysis to Encode Project
Aim 1: Haib39bioChain
- Aim 1: methylation haplotype block calling with RRBS dataset
Background
- 1, Encode: http://genome.ucsc.edu/ENCODE/downloads.html
- 2, Encode|RRBS1: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/
- 3, Encode|RRBS2: http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=437674359_aUhx08DjchWwtBjyCYv61EB7Yy8S&c=chr1&g=wgEncodeHaibMethylRrbs
- big problem for the fastq files from the above link. 1) not fastq but rather RAS file 2) reads is not true reads.
- 4, Encode|Methy450K: http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethyl450/
- 5, Encode|software: https://www.encodeproject.org/software/
Method
- 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
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
Aim 2: All the RRBS Data
- Download file
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/files.txt
- md5sum file
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibMethylRrbs/md5sum.txt
- Download fastq file with perl script
perl fastqDownloadRRBSEncode.pl files.txt
- Fastq to BAM
- BAM to hapinfo
- 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
Aim 3: All the RRBS Data (bed11)
- bed11 to bedGraph
for i in `ls *bed.gz` do zcat $i | grep -v "^track" |sort -k1,1 -k2,2n | awk '$5>9 {print $1"\t"$2"\t"$3"\t"$11}'> $i.bedGraph done
- bedGraph to bigwig
Data Download
Fastq and bed files can be downloaded from Encode Project. 101 RRBS data (bed files) were downloaded. 2,646,999 CpG loci were covered by 101 RRBS data while 866,979 CpG loci (32.8%) were detected in at least 80% samples.