Noi/NOTES/2014-5-24: Difference between revisions

From ZhangLabWiki
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'''


==== Directional mapping ====
 
* 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


==== 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 removed table because I paste the number to the wrong samples.
* I need to double check if this is the right way to get these numbers. I also need to learn how to get these files into bed file and calculate coverage at different read depth
* 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
* One thing I notice from above table is that numbers of unique CpG (1X) are similar within the same pool as shown in the same color. I am not quite sure if this is the case.
* 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 reads
apr10_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