Matt:LabNotes/2013-8-21
Jump to navigation
Jump to search
Agi26k Capture Analysis (MiSeq_130729)[edit]
- Analysis here is mostly similar to CA12k Capture (MiSeq_130325) Analysis
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 |