Matt:LabNotes/2013-8-20

From ZhangLabWiki
Jump to navigation Jump to search

CA12k Capture Analysis (MiSeq_130325)[edit]

  • Here is improved analysis of the MiSeq_130325 data
    • Previously used Bowtie2 custom parameters, which now I know are not optimal
    • Previously aligned to hg19 and now aligned to probelist used to order the probes
    • The old analysis can be found in the commented out section at the bottom of the "edit" page or Matt:LabNotes/MiSeq_130325_Analysis

Building bowtie2 index for CA12k probes[edit]

  • First converted probelist to fasta file with H1, H2,and 100bp linking region rearranged using CAProbelist_to_fasta.pl: Feb2013_200bp_probes_to_order.txt -> CAprobes_to_order_H1H2.fa
    • Fasta reference made to match reads: RevComp(H1) + RevComp(H2) + RevComp(Barcode)
  • Build index: genome-miner:$ bowtie2-build CAprobes_to_order_H1H2.fa CAprobes_H1H2

Processing reads[edit]

  • Shorten Read1 reads to 146bp to match reference length with h1h2_first146_only.pl (nothing fancy, just substr($_,0,146) for base call and base quality strings)
 D1-cDNA-RNaseA_S2_L001_R1_001.fastq -> cDNA-RNaseA_R1_146shortened.fastq
 D1-cDNAwRNaseA_S3_L001_R1_001.fastq -> cDNAwRNaseA_R1_146shortened.fastq
 D1-gDNA_S1_L001_R1_001.fastq -> gDNA_R1_146shortened.fastq

Map reads with default Bowtie2[edit]

 bowtie2 --phred33 -x CAprobes_H1H2 -q gDNA_R1_146shortened.fastq > gDNA_R1_146shortened_CAprobesH1H2.sam &
 4476041 reads; of these:
 4476041 (100.00%) were unpaired; of these:
   129983 (2.90%) aligned 0 times
   4322633 (96.57%) aligned exactly 1 time
   23425 (0.52%) aligned >1 times
 97.10% overall alignment rate
 bowtie2 --phred33 -x CAprobes_H1H2 -q cDNA-RNaseA_R1_146shortened.fastq > cDNA-RNaseA_R1_146shortened_CAprobesH1H2.sam 2> cDNA-RNase_stderr.txt &
 2852904 reads; of these:
 2852904 (100.00%) were unpaired; of these:
   81850 (2.87%) aligned 0 times
   2740392 (96.06%) aligned exactly 1 time
   30662 (1.07%) aligned >1 times
 97.13% overall alignment rate
 bowtie2 --phred33 -x CAprobes_H1H2 -q cDNAwRNaseA_R1_146shortened.fastq > cDNAwRNaseA_R1_146shortened_CAprobesH1H2.sam 2> cDNAwRNase_stderr.txt &
 2091417 reads; of these:
 2091417 (100.00%) were unpaired; of these:
   55021 (2.63%) aligned 0 times
   2010743 (96.14%) aligned exactly 1 time
   25653 (1.23%) aligned >1 times
 97.37% overall alignment rate

Sort and filter (MAPQ>=8) aligned reads[edit]

  • MAPQ=8 corresponds to 16% chance alignment is wrong, but manual checking showed that alignments with MAPQ=8 matched well
 samtools view -bS gDNA_R1_146shortened_CAprobesH1H2.sam | samtools sort - gDNA_R1_146shortened_CAprobesH1H2_sorted
 samtools view -h -F 4 -q 8 gDNA_R1_146shortened_CAprobesH1H2_sorted.bam > gDNA_R1_146shortened_CAprobesH1H2_sorted_filtered.sam
 samtools view -bS cDNA-RNaseA_R1_146shortened_CAprobesH1H2.sam | samtools sort - cDNA-RNaseA_R1_146shortened_CAprobesH1H2_sorted
 samtools view -h -F 4 -q 8 cDNA-RNaseA_R1_146shortened_CAprobesH1H2_sorted.bam > cDNA-RNaseA_R1_146shortened_CAprobesH1H2_sorted_filtered.sam
 samtools view -bS cDNAwRNaseA_R1_146shortened_CAprobesH1H2.sam | samtools sort - cDNAwRNaseA_R1_146shortened_CAprobesH1H2_sorted
 samtools view -h -F 4 -q 8 cDNAwRNaseA_R1_146shortened_CAprobesH1H2_sorted.bam > cDNAwRNaseA_R1_146shortened_CAprobesH1H2_sorted_filtered.sam

Count reads aligned for each probe[edit]

  • Used Count_reads_from_bowtie2_SAM.pl to count number of appearances of each probe in column 3 (Reference Name) of sam file
    • CountofcDNA-RNaseA_CAprobesH1H2_filtered.txt -> total count for all probes: 2,554,856
    • CountofcDNAwRNaseA_CAprobesH1H2_filtered.txt -> total count for all probes: 1,881,008
    • CountofgDNA_CAprobesH1H2_filtered.txt -> total count for all probes: 3,976,822

Histograms and Results[edit]

  • Only used cDNA-RNaseA since previous regression analysis showed with or without RNaseA were essentially the same
  • "Good" probes were the set of probes that had a capture efficiency within one order of magnitude
  • "Brain optimized" probes were the set of probes that captured

File:CA12k gDNA Capture.JPGFile:CA12k cDNA-RNaseA Capture.JPG

  • 3,141 Unique "good" genes
  • 2,615 Unique "brain optimized" genes
  • 4,185 "Brain optimized" and "good" probes
    • 1,963 "Brain optimized" and "good" genes
    • 1,849 of which are brain genes
    • 189 of which are stem cell genes

Details of genes/probes can be found in CA12k_Capture_Analysis.xlsx

Comparison of CA12k to Agi26k Results Summary[edit]

' CA12k Agi26k0gap Agi26k20gap
Good probes 7,900 6,596 5,403
Good genes 3,141 2,628 2,529
Brain-optimized probes 6,412 8,640 8,438
Brain-optimized genes 2,615 3,140 3,133
Good & Brain-optimized probes 4,185 4,822 4,006
Good & Brain-optimized genes 1,963 2,075 2,075
Good & Brain-optimized brain genes 1,849 2,002 1,996
Good & Brain-optimized stem cell genes 189 138 151