Noi/NOTES/2014-7-20

From ZhangLabWiki
Jump to navigation Jump to search

Mapping of MspI/TaqI double-digested RRBS libraries[edit]

Calendar

  • Library prep
  • Miseq Run in Palsson's lab.
    • Read length PE 150+6+150 bp
    • Data was transferred to genome-miner: /media/LTS_33T/SeqStore2/140630_MiSeq
    • Phred score: 33
  • There was a delay of getting data of this run because there were illegal characters in Palsson's lab sample IDs.
  • After I did mapping of RRBS data in previous batch, it looked like using longer adaptor sequence for trimming by Trimgalore gave higher correlation. Having discussion with Dinh, she said that we should use the longer version for trimming RRBS data generated from Illumina TruSeq adaptor with barcodes.
Read1: 
Trim read 1 with A + the “TruSeq Indexed Adapter”. Use only the prefix of the adapter sequence that is common to all Indexed Adapter sequences:
AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC

Read2:
Trim read 2 with the reverse complement of the ”TruSeq Universal Adapter”:
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT
  • The recognition site of MspI is C^CGG and TaqI is T^CGA. They both generate the same filled-in sequence (GC) at 3' end. There is no need to add extra parameter for sequence trimming since --rrbs option should deal with the filled-in sequences in both cases.

Quick data info[edit]

Sequencing IDs Libraries ID # of reads
NP-dgRRBS-6T1-Jul11-16_S1_L001_R1_001.fastq 6-T-1_MspI 5,656,056
NP-dgRRBS-6T1-Jul11-19_S2_L001_R1_001.fastq 6-T-1_MspI-TaqI 6,563,713
  • To have equal number of reads for CpG coverage comparison, I randomly picked a subset of fastq reads of index 19 using available online script subsample.py

Usage

python subsample.py <fraction> <input file 1> <input file 2> <output file 1> <output file 2>

<fraction> is a number between 0 and 1, giving the sampling faction

python subsample.py 0.855 ../NP-dgRRBS-6T1-Jul11-19_S2_L001_R1_001_val_1.fq ../NP-dgRRBS-6T1-Jul11-19_S2_L001_R2_001_val_2.fq ./sub.NP-dgRRBS-6T1-Jul11-19_S2_L001_R1_001_val_1.fq ./sub.NP-dgRRBS-6T1-Jul11-19_S2_L001_R2_001_val_2.fq &
  • I finally got 5,258,000 reads of NP-dgRRBS-6T1-Jul11-19_S2
  • # of read of sample NP-dgRRBS-6T1-Jul11-16_S1 after trimming is 5,258,311 reads
  • May need to check if this is the right way since I only subsampling one sample to have the number of read equal to the lower one.

Mapping Summary[edit]

Sample ID File ID Total reads (after COPE) Total BP # BP mapped % BP mapped CpG count (10X) CpG count (5X) CpG count (1X)
6-T-1_sRRBS NP-dgRRBS-6T1-Jul11-16_S1 9,124,296 1,078,904,486 711,246,293 65.92 1,156,962 1,692,494 2,960,140
6-T-1_dRRBS sub.NP-dgRRBS-6T1-Jul11-19_S2 9,017,527 994,374,650 584,946,934 58.83 900,920 1,557,354 3,060,748
  • Mapping rate of dRRBS is slightly lower than sRRBS. Even it's only one library from each, I think it's possible that I did not remove the band between 198-206 enough, so the repetitive sequences still existed in the library.
  • Base on the same number of read, I did not expect to get more called CpGs from dRRBS than in sRRBS using minimum read depth of 10. The reason is that dRRBS would generate more unique reads than sRRBS, so if the amount of sequencing is not high enough, they would be discarded as the read depth would be lower than 10X. To verify that

CGI annotation[edit]

  • I download CGI coordinate from genome browser
  • I then used bedtoos to find the genomic region of CpGs whether is in CGI or out side CGI. For now, I did not count CpG sites in CpG shores or CpG shelf.
Sample Total CpG (1X) # of CpG in CGI # of CpG in non-CGI % CpG in CGI % CpG in non-CGI
6-T-1_sRRBS 2,960,140 1,122,041 1,838,099 37.90 62.10
6-T-1_dRRBS 3,060,747 1,066,172 1,994,575 34.83 65.17
Sample Total CpG (10X) # of CpG in CGI # of CpG in non-CGI % CpG in CGI % CpG in non-CGI
6-T-1_sRRBS 1,156,962 498,589 658,373 43.09 56.91
6-T-1_dRRBS 900,920 379,819 521,101 42.16 57.84
  • I simply cound CpG in non-CGI region by subtracting total CpGs with CpG in CGI regions.
  • With minimum 10X the percentage of CpG in CGI from total CpGs is very close to the number reported by Lee YK, 2014.
  • Lee YK, 2014 did not report that how many reads from sRRBS and dRRBS were used to compare. They said that they can call ~1.8M CpG (min 10X) from 2 lanes of GAIIx by doing PE 36bp sequencing. They obtained ~25M reads in average or each lane. I think that they did ~50 M reads to get this data.
  • The percentage of CpG in CGI and non-CGI regions are not much different. However, as I mentioned before that number of read is not high enough to make the conclusion. I may need to look that dRRBS can gain more different CpGs compare to the data of 6-T-1 sample generated on 20140709.

# of overlapped CpGs between sRRBS (5M reads or 10M reads) and dRRBS (5M reads), minimum RD 10X[edit]

Sample Total CpGs 6-T-1_1 6-T-1_2 6-T-1_sRRBS 6-T-1_dRRBS
6-T-1_1 1,237,103 XX 1,178,553 793,314 593,142
6-T-1_2 1,562,053 1,178,553 XX 981,583 654,129
6-T-1_sRRBS 1,156,962 793,314 981,583 XX 525,240
6-T-1_dRRBS 900,920 593,142 654,129 525,240 XX