Matt:LabNotes/2013-8-21: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Mzcai
(Created page with "==Agi26k Capture Analysis (MiSeq_130729)== *Analysis here is mostly similar to CA12k Capture (MiSeq_130325) Analysis ===Building bowtie2 index for...")
 
>Mzcai
No edit summary
Line 3: Line 3:


===Building bowtie2 index for Agi26k probes===
===Building bowtie2 index for Agi26k probes===
*First converted probelist to fasta file with H1, H2,and 100bp linking region rearranged using [[Media:CAProbelist_to_fasta.txt | CAProbelist_to_fasta.pl]]: Feb2013_200bp_probes_to_order.txt -> CAprobes_to_order_H1H2.fa
*For 0gap probes converted probelist to fasta file with H1, H2,and barcode region rearranged similar to CA12kprobes [[Media:Agi26k0gapProbelist_to_fasta.txt | 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)
**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
**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 [[Media:Agi26k20gapProbelist_to_fasta_corrected.txt | 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


===Processing reads===
*Build index:
*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)
  genome-miner:$ bowtie2-build Agi26k0gapprobes_to_order_H1H2.fa Agi26k0gap_H1H2
  D1-cDNA-RNaseA_S2_L001_R1_001.fastq -> cDNA-RNaseA_R1_146shortened.fastq
  genome-miner:$ bowtie2-build Agi26k20gapprobes_to_order_H1H2_corrected.fa Agi26k20gap_H1targetH2
  D1-cDNAwRNaseA_S3_L001_R1_001.fastq -> cDNAwRNaseA_R1_146shortened.fastq
 
  D1-gDNA_S1_L001_R1_001.fastq -> gDNA_R1_146shortened.fastq
===Reads===
MZC-Agi26k-0cDNA-Jul11_CCTATT_L001_R1_001.fastq<br>
MZC-Agi26k-0gDNA-Jul11_ACTACG_L001_R1_001.fastq<br>
MZC-Agi26k-20cDNA-Jul11_CGTGTT_L001_R1_001.fastq<br>
MZC-Agi26k-20gDNA-Jul11_ACGGGT_L001_R1_001.fastq<br>


===Map reads with default Bowtie2===
===Map reads with default Bowtie2===
   bowtie2 --phred33 -x CAprobes_H1H2 -q gDNA_R1_146shortened.fastq > gDNA_R1_146shortened_CAprobesH1H2.sam &
   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 &
   4476041 reads; of these:
   4535752 reads; of these:
   4476041 (100.00%) were unpaired; of these:
   4535752 (100.00%) were unpaired; of these:
     129983 (2.90%) aligned 0 times
     81030 (1.79%) aligned 0 times
     4322633 (96.57%) aligned exactly 1 time
     4433426 (97.74%) aligned exactly 1 time
     23425 (0.52%) aligned >1 times
     21296 (0.47%) aligned >1 times
   97.10% overall alignment rate
   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 CAprobes_H1H2 -q cDNA-RNaseA_R1_146shortened.fastq > cDNA-RNaseA_R1_146shortened_CAprobesH1H2.sam 2> cDNA-RNase_stderr.txt &
   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 &
   2852904 reads; of these:
   3789673 reads; of these:
   2852904 (100.00%) were unpaired; of these:
   3789673 (100.00%) were unpaired; of these:
     81850 (2.87%) aligned 0 times
     948660 (25.03%) aligned 0 times
     2740392 (96.06%) aligned exactly 1 time
     2819628 (74.40%) aligned exactly 1 time
     30662 (1.07%) aligned >1 times
     21385 (0.56%) aligned >1 times
   97.13% overall alignment rate
   74.97% overall alignment rate


   bowtie2 --phred33 -x CAprobes_H1H2 -q cDNAwRNaseA_R1_146shortened.fastq > cDNAwRNaseA_R1_146shortened_CAprobesH1H2.sam 2> cDNAwRNase_stderr.txt &
   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 &
   2091417 reads; of these:
   5808711 reads; of these:
   2091417 (100.00%) were unpaired; of these:
   5808711 (100.00%) were unpaired; of these:
     55021 (2.63%) aligned 0 times
     873070 (15.03%) aligned 0 times
     2010743 (96.14%) aligned exactly 1 time
     4880686 (84.02%) aligned exactly 1 time
     25653 (1.23%) aligned >1 times
     54955 (0.95%) aligned >1 times
   97.37% overall alignment rate
   84.97% overall alignment rate


===Sort and filter (MAPQ>=8) aligned reads===
===Sort and filter (MAPQ>=8) aligned reads===
*MAPQ=8 corresponds to 16% chance alignment is wrong, but manual checking showed that alignments with MAPQ=8 matched well
*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 -bS Agi26k-0gDNA_R1_H1H2.sam | samtools sort - Agi26k-0gDNA_R1_H1H2_sorted
   samtools view -h -F 4 -q 8 gDNA_R1_146shortened_CAprobesH1H2_sorted.bam > gDNA_R1_146shortened_CAprobesH1H2_sorted_filtered.sam
   samtools view -h -F 4 -q 8 Agi26k-0gDNA_R1_H1H2_sorted.bam > Agi26k-0gDNA_R1_H1H2_sorted_filtered.sam


   samtools view -bS cDNA-RNaseA_R1_146shortened_CAprobesH1H2.sam | samtools sort - cDNA-RNaseA_R1_146shortened_CAprobesH1H2_sorted
   samtools view -bS Agi26k-0cDNA_R1_H1H2.sam | samtools sort - Agi26k-0cDNA_R1_H1H2_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 -h -F 4 -q 8 Agi26k-0cDNA_R1_H1H2_sorted.bam > Agi26k-0cDNA_R1_H1H2_sorted_filtered.sam


   samtools view -bS cDNAwRNaseA_R1_146shortened_CAprobesH1H2.sam | samtools sort - cDNAwRNaseA_R1_146shortened_CAprobesH1H2_sorted
   samtools view -bS Agi26k-20gDNA_R1_H1targetH2.sam | samtools sort - Agi26k-20gDNA_R1_H1targetH2_sorted
   samtools view -h -F 4 -q 8 cDNAwRNaseA_R1_146shortened_CAprobesH1H2_sorted.bam > cDNAwRNaseA_R1_146shortened_CAprobesH1H2_sorted_filtered.sam
   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===
===Count reads aligned for each probe===
*Used [[Media:Count_reads_from_bowtie2_SAM.txt | Count_reads_from_bowtie2_SAM.pl]] to count number of appearances of each probe in column 3 (Reference Name) of sam file
*Used [[Media:Count_reads_from_bowtie2_SAM.txt | Count_reads_from_bowtie2_SAM.pl]] to count number of appearances of each probe in column 3 (Reference Name) of sam file

Revision as of 23:03, 21 August 2013

Agi26k Capture Analysis (MiSeq_130729)

Building bowtie2 index for Agi26k probes

  • 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

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

 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

  • 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

  • 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

  • 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