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

From ZhangLabWiki
Jump to navigation Jump to search
>Noi
>Noi
mNo edit summary
Line 62: Line 62:
  '''Some examples of fastQC images before and after quality trimming'''
  '''Some examples of fastQC images before and after quality trimming'''
   
   
  '''10 nuclei #1'''
'''1nu #1, <u>data from lane 1</u>'''
per_base_quality, per_sequence_quality
'''Before trimming'''
[[File:apr10_1to50_no1_per_base_quality.png| 450px]] [[File:apr10_1to50_no1_per_sequence_quality.png| 450px]]
'''After trimming'''
[[File:apr10_1to50_no1_per_base_quality_trimmed.png| 450px]] [[File:apr10_1to50_no1_per_sequence_quality_trimmed.png| 450px]]
  '''10 nuclei #1, <u>data from lane 2</u>'''
  per_base_quality, per_sequence_quality
  per_base_quality, per_sequence_quality
   
   

Revision as of 23:43, 25 May 2014

RRBS data analysis

  • Data information
    • scRRBS and 10nuclei RRBS data generated on 2014-04-25 [[1]]
    • Working directory on genome-miner: /home/nplongth/Noi_scratch/2014_works/RRBS_data/RRBS_140514_HL166

Step 1: Adapter trimming and quality filtering

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

  • 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 stringency!!)
- 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

  • 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

Tool: Bismark

Step 2.1: Running bismark_genome_preparation

  • 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

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)

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
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
/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 1170863 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

  • 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
Samples Unique CpG (1x)
apr10_1to50_no1 24,776
apr10_1to50_no2 20,057
apr10_1to50_no3 24,694
apr10_1to50_no4 28,931
apr10_1to50_no5 329,622
apr10_1to50_no6 381,379
apr10_1to100_no7 78,486
apr10_1to100_no8 80,656
apr10_1to100_no9 94,957
apr10_1to100_no10 95,734
10nu_apr10_1to50_no1 118,478
10nu_apr10_1to50_no2 73,475
  • 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
  • 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.
    • 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.

Data from Tagng's group

  • 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.
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