Matt:LabNotes/2013-8-21

From ZhangLabWiki
Jump to navigation Jump to search

Agi26k Capture Analysis (MiSeq_130729)[edit]

Building bowtie2 index for Agi26k probes[edit]

  • For 0gap probes converted probelist to fasta file with H1, H2,and barcode region rearranged similar to CA12kprobes Agi26k0gapProbelist_to_fasta.pl: FISSEQ_Probes_Dec2012.txt -> Agi26k0gapprobes_to_order_H1H2.fa
    • Fasta reference made to match reads: RevComp(H1) + RevComp(H2) + RevComp(Barcode)
    • Had to delete probes 12964-26142 since they are 20gap probes
  • For 20gap probes did something similar but had to include 20bp target region in between H1 and H2 Agi26k20gapProbelist_to_fasta_corrected.pl: FISSEQ_Probes_Dec2012.txt -> Agi26k20gapprobes_to_order_H1H2_corrected.fa
    • Target sequence was found by using H1 sequence to search through outputFile_20gap.txt and outputFile_20gap_stem.txt (the two files output from ppDesigner when designing probes for genes and stem cell genes)
    • Fasta reference made to match reads: RevComp(H1) + target + RevComp(H2) + RevComp(Barcode)
    • Had to delete probes 0-12963 since they are 0gap probes
  • Build index:
 genome-miner:$ bowtie2-build Agi26k0gapprobes_to_order_H1H2.fa Agi26k0gap_H1H2
 genome-miner:$ bowtie2-build Agi26k20gapprobes_to_order_H1H2_corrected.fa Agi26k20gap_H1targetH2

Reads[edit]

MZC-Agi26k-0cDNA-Jul11_CCTATT_L001_R1_001.fastq
MZC-Agi26k-0gDNA-Jul11_ACTACG_L001_R1_001.fastq
MZC-Agi26k-20cDNA-Jul11_CGTGTT_L001_R1_001.fastq
MZC-Agi26k-20gDNA-Jul11_ACGGGT_L001_R1_001.fastq

Map reads with default Bowtie2[edit]

 bowtie2 --phred33 -x Agi26k0gap_H1H2 -q MZC-Agi26k-0gDNA-Jul11_ACTACG_L001_R1_001.fastq > Agi26k-0gDNA_R1_H1H2.sam 2> Agi26k-0gDNA_stderr.txt &
 4535752 reads; of these:
 4535752 (100.00%) were unpaired; of these:
   81030 (1.79%) aligned 0 times
   4433426 (97.74%) aligned exactly 1 time
   21296 (0.47%) aligned >1 times
 98.21% overall alignment rate
 bowtie2 --phred33 -x Agi26k0gap_H1H2 -q MZC-Agi26k-0cDNA-Jul11_CCTATT_L001_R1_001.fastq > Agi26k-0cDNA_R1_H1H2.sam 2> Agi26k-0cDNA_stderr.txt &
 4981156 reads; of these:
 4981156 (100.00%) were unpaired; of these:
   84336 (1.69%) aligned 0 times
   4836874 (97.10%) aligned exactly 1 time
   59946 (1.20%) aligned >1 times
 98.31% overall alignment rate
 bowtie2 --phred33 -x Agi26k20gap_H1targetH2 -q MZC-Agi26k-20gDNA-Jul11_ACGGGT_L001_R1_001.fastq > Agi26k-20gDNA_R1_H1targetH2.sam 2> Agi26k-20gDNA_stderr.txt &
 3789673 reads; of these:
 3789673 (100.00%) were unpaired; of these:
   948660 (25.03%) aligned 0 times
   2819628 (74.40%) aligned exactly 1 time
   21385 (0.56%) aligned >1 times
 74.97% overall alignment rate
 bowtie2 --phred33 -x Agi26k20gap_H1targetH2 -q MZC-Agi26k-20cDNA-Jul11_CGTGTT_L001_R1_001.fastq > Agi26k-20cDNA_R1_H1targetH2.sam 2> Agi26k-20cDNA_stderr.txt &
 5808711 reads; of these:
 5808711 (100.00%) were unpaired; of these:
   873070 (15.03%) aligned 0 times
   4880686 (84.02%) aligned exactly 1 time
   54955 (0.95%) aligned >1 times
 84.97% 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 Agi26k-0gDNA_R1_H1H2.sam | samtools sort - Agi26k-0gDNA_R1_H1H2_sorted
 samtools view -h -F 4 -q 8 Agi26k-0gDNA_R1_H1H2_sorted.bam > Agi26k-0gDNA_R1_H1H2_sorted_filtered.sam
 samtools view -bS Agi26k-0cDNA_R1_H1H2.sam | samtools sort - Agi26k-0cDNA_R1_H1H2_sorted
 samtools view -h -F 4 -q 8 Agi26k-0cDNA_R1_H1H2_sorted.bam > Agi26k-0cDNA_R1_H1H2_sorted_filtered.sam
 samtools view -bS Agi26k-20gDNA_R1_H1targetH2.sam | samtools sort - Agi26k-20gDNA_R1_H1targetH2_sorted
 samtools view -h -F 4 -q 8 Agi26k-20gDNA_R1_H1targetH2_sorted.bam > Agi26k-20gDNA_R1_H1targetH2_sorted_filtered.sam
 samtools view -bS Agi26k-20cDNA_R1_H1targetH2.sam | samtools sort - Agi26k-20cDNA_R1_H1targetH2_sorted
 samtools view -h -F 4 -q 8 Agi26k-20cDNA_R1_H1targetH2_sorted.bam > Agi26k-20cDNA_R1_H1targetH2_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
    • CountofgDNA_Agi26k0gapH1H2_filtered.txt -> total count for all probes: 4,414,961
    • CountofcDNA_Agi26k0gapH1H2_filtered.txt -> total count for all probes: 4,850,618
    • CountofgDNA_Agi26k20gapH1targetH2_filtered.txt -> total count for all probes: 1,710,730
    • CountofcDNA_Agi26k20gapH1targetH2_filtered.txt -> total count for all probes: 4,298,273

Histograms and Results[edit]

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

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

Agi26k0gap[edit]

File:Agi26k0gap gDNA Capture.JPGFile:Agi26k0gap cDNA Capture.JPG

  • 2,628 Unique "good" genes
  • 3,140 Unique "brain optimized" genes
  • 4,822 "Brain optimized" and "good" probes
    • 2,075 "Brain optimized" and "good" genes
    • 2,002 of which are brain genes
    • 138 of which are stem cell genes

Agi26k20gap[edit]

File:Agi26k20gap gDNA Capture.JPGFile:Agi26k20gap cDNA Capture.JPG

  • 2,529 Unique "good" genes
  • 3,133 Unique "brain optimized" genes
  • 4,006 "Brain optimized" and "good" probes
    • 2,075 "Brain optimized" and "good" genes
    • 1,996 of which are brain genes
    • 151 of which are stem cell genes

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