Noi/NOTES/2014-5-24: Difference between revisions
Jump to navigation
Jump to search
>Noi mNo edit summary |
>Noi |
||
(4 intermediate revisions by the same user not shown) | |||
Line 1: | Line 1: | ||
= RRBS data analysis = | = RRBS data analysis = | ||
'''Link to calendar: [[http://genome-tech.ucsd.edu/LabNotes/index.php/noi:DMR220k_LabNotes#2014]]''' | |||
* Data information | * Data information | ||
** scRRBS and 10nuclei RRBS data generated on 2014-04-25 [[http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2014-4-25]] | ** scRRBS and 10nuclei RRBS data generated on 2014-04-25 [[http://genome-tech.ucsd.edu/LabNotes/index.php/Noi/NOTES/2014-4-25]] | ||
Line 129: | Line 130: | ||
* It generated '''Bisulfite_Genome''' folder, including two sub-folders '''CT_conversion''' and '''GA_conversion''' | * It generated '''Bisulfite_Genome''' folder, including two sub-folders '''CT_conversion''' and '''GA_conversion''' | ||
==== Step 2.2: Running bismark mapping (using bowtie) ==== | ==== Step 2.2: Running bismark mapping (using bowtie) ==== | ||
==== Directional mapping, bowtie ==== | |||
* In folder '''NP.RRBS.Trim1_directional/Bismark-run.directional''' | |||
'''Run.bismark.directional.sh''' | '''Run.bismark.directional.sh''' | ||
Line 267: | Line 260: | ||
* I simply use awk command to get chromosome position and start position in column $3 and $4, respectively in CpG_context.XXXX._bismark.txt files. Then sort out the unique line to get minimum 1x unique CpGs | * I simply use awk command to get chromosome position and start position in column $3 and $4, respectively in CpG_context.XXXX._bismark.txt files. Then sort out the unique line to get minimum 1x unique CpGs | ||
Ex. awk '{print $3"\t"$4}' CpG_context_apr10_1to50_no1.txt_trimmed.fq_bismark.txt | awk '{if($1~/chr/) print $0}' | sort -u | wc -l | Ex. awk '{print $3"\t"$4}' CpG_context_apr10_1to50_no1.txt_trimmed.fq_bismark.txt | awk '{if($1~/chr/) print $0}' | sort -u | wc -l | ||
* I also need to run with bowtie2 to compare the result with bowtie1 | * I also need to run with bowtie2 to compare the result with bowtie1 | ||
* Need to double check and learn bismark in great details. | * Need to double check and learn bismark in great details. | ||
* There is another point I need to add to libraries in pool1 (green highlighted). When I did PAGE size-selection, the gel was badly stuck to the cassette, and I may not cut at the right size and I lost some lane. The yield of this pool was very low. I am not sure if this is the reason that I got unique CpG lower than other pools. | * There is another point I need to add to libraries in pool1 (green highlighted). When I did PAGE size-selection, the gel was badly stuck to the cassette, and I may not cut at the right size and I lost some lane. The yield of this pool was very low. I am not sure if this is the reason that I got unique CpG lower than other pools. | ||
Line 285: | Line 275: | ||
:- I have not continued to extract methylation data into bed file and find out the coverage at different read depth yet. I need your suggestion before go to the next step. | :- I have not continued to extract methylation data into bed file and find out the coverage at different read depth yet. I need your suggestion before go to the next step. | ||
:''Please go ahead to extract the methylation information, and share the bed files with me. I’d like to take a look by myself in next day or two.'' | :''Please go ahead to extract the methylation information, and share the bed files with me. I’d like to take a look by myself in next day or two.'' | ||
=== Summary of CpG coverage at different read dept === | === Summary of CpG coverage at different read dept === | ||
* After running bismark_methylation_extractor, I generate bed file from CpG_context out file | * After running bismark_methylation_extractor, I generate bed file from CpG_context out file | ||
/home/nplongth/softwares/bismark_v0.12.2/bismark2bedGraph $f -o $f.bed | /home/nplongth/softwares/bismark_v0.12.2/bismark2bedGraph $f -o $f.bed | ||
<u>input file:</u> $f = output file from bismark_methylation_extractor | <u>input file:</u> $f = output file from bismark_methylation_extractor | ||
* Header of bed file: <chromosome> <start position> <end position> <methylation percentage> <count methylated> <count unmethylated> | <u>output file:</u> $f.bed & $f.bed.bismark.cov (includes 2 more columns of methylated and unmethylated count) | ||
* Header of bed file (original file is $f.cov file) | |||
: <chromosome> <start position> <end position> <methylation percentage> <count methylated> <count unmethylated> | |||
{| {{table}} class =wikitable | {| {{table}} class =wikitable | ||
| align="center" style="width:120px; background:#f0f0f0;"|'''Samples''' | | align="center" style="width:120px; background:#f0f0f0;"|'''Samples''' |
Latest revision as of 04:00, 3 June 2014
RRBS data analysis[edit]
Link to calendar: [[1]]
- Data information
- scRRBS and 10nuclei RRBS data generated on 2014-04-25 [[2]]
- Working directory on genome-miner: /home/nplongth/Noi_scratch/2014_works/RRBS_data/RRBS_140514_HL166
Step 1: Adapter trimming and quality filtering[edit]
Tool: Trim Galore from Babrahan Bioinformatic which is a wrapper script that use Cutadapt for adapter trimming and fastQC for quality filtering
- By Running checkphred.pl, I confirm that all sequencing data are phred64
Directional option[edit]
- In folder NP.RRBS.Trim1_directional
- I am not authorized to install python-pip, so there were some difficulties to install functional cutadapt since the library was missing. My friend helped me to generate virtual environment, and it needs to be activated when running cutadapt and trim_galore.
- Activate environment: source ~/env/bin/activate
- Deactivate environment: deactivate
#!/bin/bash for f in ../RRBS_140514_HL166.rawdata/*.txt do ../trim_galore -q 20 --phred64 -a AGATCGGAAGAGC --rrbs $f done
- - Quality Phred score cutoff: 20
- - Quality encoding type selected: ASCII+64
- - Adapter sequence: 'AGATCGGAAGAGC'
- - Maximum trimming error rate: 0.1 (default)
- - Minimum required adapter overlap (stringency): 1 bp (most stringent!!)
- - Minimum required sequence length before a sequence gets removed: 20 bp
- - File was specified to be an MspI-digested RRBS sample. Sequences with adapter contamination will be trimmed a further 2 bp to remove potential methylation-biased bases from the end-repair reaction
Samples # of reads before trimming # of reads after trimming # of remaining reads apr10_1to50_no1 5,659,042 3,996,745 70.63 apr10_1to50_no2 5,707,396 3,863,720 67.70 apr10_1to50_no3 5,492,909 3,922,931 71.42 apr10_1to50_no4 5,351,954 3,526,832 65.90 apr10_1to50_no5 3,605,391 2,454,517 68.08 apr10_1to50_no6 4,027,892 2,651,359 65.82 apr10_1to100_no7 3,352,177 2,250,422 67.13 apr10_1to100_no8 2,071,082 1,428,535 68.98 apr10_1to100_no9 2,752,431 1,929,729 70.11 apr10_1to100_no10 2,203,299 1,424,776 64.67 10nu_apr10_1to50_no1 4,527,606 3,163,200 69.86 10nu_apr10_1to50_no2 4,601,068 3,012,729 65.48
File:Trimgalore directional readcompare.png Some examples of fastQC images before and after quality trimming 1nu #1, data from lane 1 per_base_quality, per_sequence_quality Before trimming File:Apr10 1to50 no1 per base quality.png File:Apr10 1to50 no1 per sequence quality.png After trimming File:Apr10 1to50 no1 per base quality trimmed.png File:Apr10 1to50 no1 per sequence quality trimmed.png 10 nuclei #1, data from lane 2 per_base_quality, per_sequence_quality Before trimming File:10nu apr10 1to50 no1 per base quality.png File:10nu apr10 1to50 no1 per sequence quality.png After trimming File:10nu apr10 1to50 no1 per base quality trimmed.png File:10nu apr10 1to50 no1 per sequence quality trimmed.png
Non-directional option[edit]
- I also ran trim_galore with --non_directional option. "This option will screen quality-trimmed sequences for 'CAA' or 'CGA' at the start of the read and, if found, removes the first two basepairs. Like with the option '--rrbs' this avoids using cytosine positions that were filled-in during the end-repair step. '--non_directional' requires '--rrbs' to be specified as well" (Krueger F, 2011).
- In folder NP.RRBS.Trim2_non-directional
#!/bin/bash for f in ../RRBS_140514_HL166.rawdata/*.txt do ../trim_galore -q 20 --phred64 -a AGATCGGAAGAGC --rrbs --non_directional $f done
Samples # of reads before trimming # of reads after trimming # of remaining reads apr10_1to50_no1 5,659,042 3,996,732 70.63 apr10_1to50_no2 5,707,396 3,863,706 67.70 apr10_1to50_no3 5,492,909 3,922,916 71.42 apr10_1to50_no4 5,351,954 3,526,823 65.90 apr10_1to50_no5 3,605,391 2,454,507 68.08 apr10_1to50_no6 4,027,892 2,651,348 65.82 apr10_1to100_no7 3,352,177 2,250,419 67.13 apr10_1to100_no8 2,071,082 1,428,534 68.98 apr10_1to100_no9 2,752,431 1,929,725 70.11 apr10_1to100_no10 2,203,299 1,424,771 64.67 10nu_apr10_1to50_no1 4,527,606 3,163,194 69.86 10nu_apr10_1to50_no2 4,601,068 3,012,718 65.48
- I got exact the same number of remaining reads as directional option
Step 2: Mapping[edit]
Tool: Bismark
Step 2.1: Running bismark_genome_preparation[edit]
- In directory: /home/nplongth/hg19_Bismark/, I copied hg19 reference genome sequences (HsGenome19.fa) from /home/kunzhang/HsGenome/hg19/
../softwares/bismark_v0.12.2/bismark_genome_preparation --path_to_bowtie /home/kunzhang/softwares/bowtie-0.12.8/ --verbose /home/nplongth/hg19_Bismark/
- It generated Bisulfite_Genome folder, including two sub-folders CT_conversion and GA_conversion
Step 2.2: Running bismark mapping (using bowtie)[edit]
Directional mapping, bowtie[edit]
- In folder NP.RRBS.Trim1_directional/Bismark-run.directional
Run.bismark.directional.sh #!/bin/bash for f in ../*.fq do /home/nplongth/softwares/bismark_v0.12.2/bismark --path_to_bowtie /home/nplongth/softwares/bowtie-0.12.8 -n 1 /home/nplongth/hg19_Bismark/ $f done
Extract number from bismark report files
Bowtie, directional option[edit]
Samples # of raw reads # of read analyzed
# of read with unique alignment Mapping efficiency (%) Mapping efficiency
to raw readsapr10_1to50_no1 5,659,042 3,996,745 2,371,651 59.34 41.91 apr10_1to50_no2 5,707,396 3,863,720 2,129,134 55.11 37.30 apr10_1to50_no3 5,492,909 3,922,931 2,024,895 51.62 36.86 apr10_1to50_no4 5,351,954 3,526,832 1,830,561 51.90 34.20 apr10_1to50_no5 3,605,391 2,454,517 1,175,117 47.88 32.59 apr10_1to50_no6 4,027,892 2,651,359 1,263,436 47.65 31.37 apr10_1to100_no7 3,352,177 2,250,422 1,268,116 56.35 37.83 apr10_1to100_no8 2,071,082 1,428,535 873,250 61.13 42.16 apr10_1to100_no9 2,752,431 1,929,729 785,467 40.70 28.54 apr10_1to100_no10 2,203,299 1,424,776 791,273 55.54 35.91 10nu_apr10_1to50_no1 4,527,606 3,163,200 1,383,758 43.75 30.56 10nu_apr10_1to50_no2 4,601,068 3,012,729 1,659,102 55.07 36.06
- - # of raw reads = # of read after trimming
- - # of read analyzed = # of remaning reads after trimming
- - I am not sure which numbers should be used to present actual mapping rate or mapping efficiency between comparing to the original raw reads and to the reads used for mapping.
Samples Ave. mapping efficiency STDEV Ave. mapping efficiency to raw reads STDEV 1 nuclei 52.72 6.12 35.87 4.35 10nuclei 49.41 8.01 33.31 3.89
Bowtie, non-directional option[edit]
/home/nplongth/softwares/bismark_v0.12.2/bismark --non_directional --path_to_bowtie /home/nplongth/softwares/bowtie-0.12.8 -n 1 /home/nplongth/hg19_Bismark/ $f
Samples # of raw reads # of read analyzed # of read with unique alignment Mapping efficiency (%) Mapping efficiency to raw reads apr10_1to50_no1 5,659,042 3,996,732 2,186,983 54.72 38.65 apr10_1to50_no2 5,707,396 3,863,706 2,005,117 51.90 35.13 apr10_1to50_no3 5,492,909 3,922,916 1,948,269 49.66 35.47 apr10_1to50_no4 5,351,954 3,526,823 1,719,546 48.76 32.13 apr10_1to50_no5 3,605,391 2,454,507 1,120,980 45.67 31.09 apr10_1to50_no6 4,027,892 2,651,348 1,214,334 45.80 30.15 apr10_1to100_no7 3,352,177 2,250,419 1,170,863 52.03 34.93 apr10_1to100_no8 2,071,082 1,428,534 826,575 57.86 39.91 apr10_1to100_no9 2,752,431 1,929,725 747,001 38.71 27.14 apr10_1to100_no10 2,203,299 1,424,771 716,387 50.28 32.51 10nu_apr10_1to50_no1 4,527,606 3,163,194 1,254,365 39.66 27.70 10nu_apr10_1to50_no2 4,601,068 3,012,718 1,475,133 48.96 32.06
Samples Ave. mapping efficiency STDEV Ave. mapping efficiency to raw reads STDEV 1 nuclei 49.54 5.33 33.71 3.89 10nuclei 44.31 6.58 29.88 3.08
- I am working on RRBS data from Tang's group to compare if I could get similar numbers by using the method I use here and the number reported in his paper.
- I then will continue to extracting CpG methylation data
Step 2.3: Running bismark methylation_extractor[edit]
- In directory: NP.RRBS.Trim1_directional/Bismark-run.directional
- I used default set up
Run.bismark.methylationExtractor.sh for f in ./*bismark.sam do /home/nplongth/softwares/bismark_v0.12.2/bismark_methylation_extractor -s --comprehensive $f done
- I simply use awk command to get chromosome position and start position in column $3 and $4, respectively in CpG_context.XXXX._bismark.txt files. Then sort out the unique line to get minimum 1x unique CpGs
Ex. awk '{print $3"\t"$4}' CpG_context_apr10_1to50_no1.txt_trimmed.fq_bismark.txt | awk '{if($1~/chr/) print $0}' | sort -u | wc -l
- I also need to run with bowtie2 to compare the result with bowtie1
- Need to double check and learn bismark in great details.
- There is another point I need to add to libraries in pool1 (green highlighted). When I did PAGE size-selection, the gel was badly stuck to the cassette, and I may not cut at the right size and I lost some lane. The yield of this pool was very low. I am not sure if this is the reason that I got unique CpG lower than other pools.
2014-05-27 : Response from Dr. Zhang[edit]
- - I used bowtie for mapping and haven't tried bowtie2. I can try it out later to see the different.
- It’s okay for now. In the future, I think we want to run adaptor trimming, then map the reads using our own bisReadMapper
- - I am not very sure if I should treat RRBS data as directional or non-directional. I have tried both ways. It turned out that non-directional option has lower mapping rate.
- I believe it’s directional, as the adaptors were added prior to bisulfite conversion.
- - For bismark mapping, I allowed 1 non-bisulfite mismatch per read. I am not sure if this is too stringent.
- Should be fine for short reads. When have longer HiSeq reads in the future, you want to allow for 1 mismatch per 30bp.
- - I have not continued to extract methylation data into bed file and find out the coverage at different read depth yet. I need your suggestion before go to the next step.
- Please go ahead to extract the methylation information, and share the bed files with me. I’d like to take a look by myself in next day or two.
Summary of CpG coverage at different read dept[edit]
- After running bismark_methylation_extractor, I generate bed file from CpG_context out file
/home/nplongth/softwares/bismark_v0.12.2/bismark2bedGraph $f -o $f.bed input file: $f = output file from bismark_methylation_extractor output file: $f.bed & $f.bed.bismark.cov (includes 2 more columns of methylated and unmethylated count)
- Header of bed file (original file is $f.cov file)
- <chromosome> <start position> <end position> <methylation percentage> <count methylated> <count unmethylated>
Samples | # of reads before trimming | # of reads after trimming | Unique CpG (1X) | Unique CpG (5X) | Unique CpG (10X) |
1nu_no1 | 5,659,042 | 3,996,732 | 118,478 | 66,644 | 53,058 |
1nu_no2 | 5,707,396 | 3,863,706 | 73,475 | 44,483 | 34,667 |
1nu_no3 | 5,492,909 | 3,922,916 | 95,734 | 51,733 | 41,433 |
1nu_no4 | 5,351,954 | 3,526,823 | 78,486 | 43,788 | 35,098 |
1nu_no5 | 3,605,391 | 2,454,507 | 80,656 | 37,718 | 29,394 |
1nu_no6 | 4,027,892 | 2,651,348 | 94,957 | 44,828 | 34,796 |
1nu_no7 | 3,352,177 | 2,250,419 | 24,776 | 13,255 | 11,172 |
1nu_no8 | 2,071,082 | 1,428,534 | 20,057 | 10,330 | 8,490 |
1nu_no9 | 2,752,431 | 1,929,725 | 24,694 | 13,425 | 11,098 |
1nu_no10 | 2,203,299 | 1,424,771 | 28,931 | 16,002 | 12,899 |
10nu_no1 | 4,527,606 | 3,163,194 | 329,622 | 160,782 | 113,091 |
10nu_no2 | 4,601,068 | 3,012,718 | 381,379 | 185,896 | 131,392 |
- From table above
- 1nu_no1- 6: I did ligation using 1ul of 1:50 diluted adaptors (blue)
- 1nu_no7-10: I did ligation using 1ul of 1:100 diluted adaptors (no color)
- 10nu_no1-2: I did ligation using 1ul of 1:50 diluted adaptors (pink)
- For 1nuclei samples, it's obvious that amount of adaptors in ligation effects the coverage of CpGs as 1:100 diluted adaptors have much lower coverage compare to 1:50 dilution. 1nu_no5,6, and 7 are a good indicator as they have pretty close sequencing amount.
- # of CpG coverage could be improved by using more concentrate adaptors in ligation reaction (at least 1:50 dilution) and increase sequencing amount and read length. The quality of sequencing is one of the factors as well.
Data from Dr. Tang's group[edit]
- Working directory on genome-miner: /home/nplongth/Noi_scratch/2014_works/RRBS_data/RRBS_140514_HL166/GuoH_scRRBSdata
- I downloaded sra files of mESC single cell #1-#8, pooled 10 cells and 20 cells, and bulk mESC and split into read1 and read2 fastq files with sratoolkit
/home/nplongth/softwares/sratoolkit.2.3.5-2-ubuntu64/bin/fastq-dump --split-files $f
- Based on data I retrieved from GEO accession #GSE47343, they did PE 101bp for the samples in the table below. I summarized only total # of raw reads and unique CpG (1x) in this table.
- I have not continued to test on this data yet. Need to wait for Dr. Zhang's opinion about the parameters I used for adapter trimming, quality filtering and mapping.
Sample # of paired-end reads Unique CpG (1x) scRRBS_mESC_Single_Cell1 15,574,909 1,309,191 scRRBS_mESC_Single_Cell2 11,592,696 955,619 scRRBS_mESC_Single_Cell3 11,072,280 1,056,351 scRRBS_mESC_Single_Cell4 19,049,260 1,535,234 scRRBS_mESC_Single_Cell5 5,183,527 1,269,763 scRRBS_mESC_Single_Cell6 9,943,188 970,525 scRRBS_mESC_Single_Cell7 5,831,225 496,715 scRRBS_mESC_Single_Cell8 6,198,312 573,049 RRBS_mESC_Pooled_10Cell 19,518,011 2,381,797 RRBS_mESC_Pooled_20Cell 19,655,365 2,592,919 RRBS_mESC_Bulk 11,240,829 2,411,401