Matt:LabNotes/2015-4-8
CA12k_Nov2014_V7 in vitro Capture Sequencing Analysis[edit]
- With in vitro Capture data, DARTFISH gene counts can be normalized and compared to other methods such as RNA-Seq
- For V4 in vitro Capture Sequencing Analysis see: Matt:LabNotes/2015-3-19
Mapping MiSeq reads to Oligo Sequences[edit]
Convert Probelist to Fasta File[edit]
CA12kNov2014_Probelist2Fasta.pl
bowtie2-build CA12k_Nov2014_V7_H1H2.fa CA12k_Nov2014_V7_H1H2
Map with Bowtie2[edit]
Mapping full read to reference[edit]
- Had a 0% mapping rate
- Reason is because the tail end of 100bp reads had very low base quality scores (some reads had up 80 "#")
- Forums call this "B-tailing" because in Phred64 'B' represents score of 2 while Phred33 uses '#'
- This means something went wrong with sequencing at a certain cycle (possibly due to fragment being shorter than read length but this library should definitely be long enough)
Use fastx toolkit to visualize read quality stats[edit]
/home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastx_quality_stats/fastx_quality_stats -Q33 -i MC-20150121_CA12kNov2014_V7gDNA-4.R1.fastq -o MC-20150121_CA12kNov2014_V7gDNA_qualstats.txt /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastx_quality_stats/fastx_quality_stats -Q33 -i MC-20150121_CA12kNov2014_V7cDNA-5.R1.fastq -o MC-20150121_CA12kNov2014_V7cDNA_qualstats.txt
/home/kunzhang/softwares/fastx_toolkit-0.0.13.2/scripts/fastq_quality_boxplot_graph.sh -i MC-20150121_CA12kNov2014_V7gDNA_qualstats.txt -o MC-20150121_CA12kNov2014_V7gDNA_qualstats.png -t CA12kNov2014_V7_gDNA /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/scripts/fastq_quality_boxplot_graph.sh -i MC-20150121_CA12kNov2014_V7cDNA_qualstats.txt -o MC-20150121_CA12kNov2014_V7cDNA_qualstats.png -t CA12kNov2014_V7_cDNA
File:MC-20150121 CA12kNov2014 V7gDNA qualstats.png File:MC-20150121 CA12kNov2014 V7cDNA qualstats.png
Trim reads with fastx toolkit[edit]
- trimmed.fastq: Trim by base quality (minimum quality=10, minimum length=20)
/home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastq_quality_trimmer/fastq_quality_trimmer -Q33 -v -t 10 -l 20 -i MC-20150121_CA12kNov2014_V7gDNA-4.R1.fastq -o MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmed.fastq Minimum Quality Threshold: 10 Minimum Length: 20 Input: 2390629 reads. Output: 2390528 reads. discarded 101 (0%) too-short reads. /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastq_quality_trimmer/fastq_quality_trimmer -Q33 -v -t 10 -l 20 -i MC-20150121_CA12kNov2014_V7cDNA-5.R1.fastq -o MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmed.fastq Minimum Quality Threshold: 10 Minimum Length: 20 Input: 3675652 reads. Output: 3675504 reads. discarded 148 (0%) too-short reads.
- trimmedv2.fastq: Trim by fixed length (keep bases 2-45; length = 44)
/home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastx_trimmer/fastx_trimmer -Q33 -f 2 -l 45 -i MC-20150121_CA12kNov2014_V7gDNA-4.R1.fastq -o MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmedv2.fastq /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastx_trimmer/fastx_trimmer -Q33 -f 2 -l 45 -i MC-20150121_CA12kNov2014_V7cDNA-5.R1.fastq -o MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmedv2.fastq
Map quality trimmed reads to reference[edit]
- Mapping rate still too low
- Looking at fastq by eye shows that the algorithm is missing easy ones that I can find by ctrl-f
bowtie2 --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmed.fastq > CA12kNov2014_V7gDNA_R1_H1H2.sam 2> CA12kNov2014_V7gDNA_stderr.txt & 2390528 reads; of these: 2390528 (100.00%) were unpaired; of these: 2199579 (92.01%) aligned 0 times 190821 (7.98%) aligned exactly 1 time 128 (0.01%) aligned >1 times 7.99% overall alignment rate bowtie2 --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmed.fastq > CA12kNov2014_V7cDNA_R1_H1H2.sam 2> CA12kNov2014_V7cDNA_stderr.txt & 3675504 reads; of these: 3675504 (100.00%) were unpaired; of these: 3366177 (91.58%) aligned 0 times 309126 (8.41%) aligned exactly 1 time 201 (0.01%) aligned >1 times 8.42% overall alignment rate
Map v2 trimmed reads to reference[edit]
- Much better mapping rate but examining fastq by eye still shows the algorithm can do better
- A common one missed is PDE1A because it has a repeating "AC" (softmasked) that I guess has a lot of sequencing error
bowtie2 --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmedv2.fastq > CA12kNov2014_V7gDNA_R1_H1H2.sam 2> CA12kNov2014_V7gDNA_stderr.txt & 2390629 reads; of these: 2390629 (100.00%) were unpaired; of these: 191725 (8.02%) aligned 0 times 2189100 (91.57%) aligned exactly 1 time 9804 (0.41%) aligned >1 times 91.98% overall alignment rate bowtie2 --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmedv2.fastq > CA12kNov2014_V7cDNA_R1_H1H2.sam 2> CA12kNov2014_V7cDNA_stderr.txt & 3675652 reads; of these: 3675652 (100.00%) were unpaired; of these: 275464 (7.49%) aligned 0 times 3385203 (92.10%) aligned exactly 1 time 14985 (0.41%) aligned >1 times 92.51% overall alignment rate
Map v2 trimmed reads with Local alignment[edit]
- Best mapping rate so far, I guess local alignment helps with all the sequencing errors in the reads?
- Overall alignment rate comparable to V4 in vitro capture alignment rate
- High multiple alignments but algorithm should choose the best mapping score
- USING THIS SAM FILE FOR FURTHER ANALYSIS
bowtie2 --local --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmedv2.fastq > CA12kNov2014_V7gDNA_R1_H1H2_localv2.sam 2> CA12kNov2014_V7gDNA_stderr_localv2.txt & 2390629 reads; of these: 2390629 (100.00%) were unpaired; of these: 76473 (3.20%) aligned 0 times 886116 (37.07%) aligned exactly 1 time 1428040 (59.73%) aligned >1 times 96.80% overall alignment rate bowtie2 --local --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmedv2.fastq > CA12kNov2014_V7cDNA_R1_H1H2_localv2.sam 2> CA12kNov2014_V7cDNA_stderr_localv2.txt & 3675652 reads; of these: 3675652 (100.00%) were unpaired; of these: 109047 (2.97%) aligned 0 times 1368908 (37.24%) aligned exactly 1 time 2197697 (59.79%) aligned >1 times 97.03% overall alignment rate
samtools view -bS CA12kNov2014_V7gDNA_R1_H1H2_localv2.sam | samtools sort - CA12kNov2014_V7gDNA_R1_H1H2_sorted samtools view -h -F 4 CA12kNov2014_V7gDNA_R1_H1H2_sorted.bam > CA12kNov2014_V7gDNA_R1_H1H2_sorted_filtered.sam
samtools view -bS CA12kNov2014_V7cDNA_R1_H1H2_localv2.sam | samtools sort - CA12kNov2014_V7cDNA_R1_H1H2_sorted samtools view -h -F 4 CA12kNov2014_V7cDNA_R1_H1H2_sorted.bam > CA12kNov2014_V7cDNA_R1_H1H2_sorted_filtered.sam
- Using local alignment for quality trimmed reads leads to too many multiple alignments
bowtie2 --local --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmed.fastq > CA12kNov2014_V7cDNA_R1_H1H2_local.sam 2> CA12kNov2014_V7cDNA_stderr_local.txt & 3675504 reads; of these: 3675504 (100.00%) were unpaired; of these: 11255 (0.31%) aligned 0 times 64099 (1.74%) aligned exactly 1 time 3600150 (97.95%) aligned >1 times 99.69% overall alignment rate
Count # of Reads Aligned to each Probe[edit]
- Output:
CountReadsPer_Gene_Probe_excludeSoftMask.pl
- Output:
- CA12kNov2014_V7_SoftMaskedProbes.txt
- All probes containing softmasked regions, '1' in 2nd column indicates it also is a part of CA12kNov2014_V4_SoftMaskedProbes.txt
- CA12kNov2014_V7gDNA_R1_H1H2_sorted_filtered_Excluded_Probecounts.txt
- Counts for only probes that are not in CA12kNov2014_V7_SoftMaskedProbes.txt
- CA12kNov2014_V7gDNA_R1_H1H2_sorted_filtered_Excluded_Genecounts.txt
- Counts for only genes that have no SoftMaskedProbes
- CA12kNov2014_V7gDNA_R1_H1H2_sorted_filtered_ExcludedProbesOnly_Genecounts
- Counts for all genes but do not include count for SoftMaskedProbes
- CA12kNov2014_V7_SoftMaskedProbes.txt
Regression Analysis[edit]
- DARTFISH Data:
- Px-px decoding of PGP1f V4 (S2 & S5)
- Px-px decoding of PGP1f V7 (S1)
- Px-px decoding of BA8 V4 (S1)
- V4 and V7 in vitro counts
- PGP1f RNA-seq Data
- PGP1f BeadArray Data
- PGP1f FISSEQ Data
Control: Variability between samples[edit]
- PGP1f V4 S2 vs PGP1f V4 S5
- Normalize counts with V4 in vitro counts and log transform
File:PGP1fV4 log norm S2 vs S5.jpeg
V4 vs V7[edit]
- PGP1f V4 (S2+S5) vs PGP1f V7 (S1)
- Log transform
File:PGP1f log sum V4 vs V7.jpeg
- Normalize counts with in vitro counts and log transform
- Low correlation is disheartening, maybe V7 in vitro counts are off
File:PGP1f log norm V4 vs V7.jpeg
PGP1f vs BA8[edit]
- PGP1f V4 (S2+S5) vs BA8 V4 (S1)
- Log transform
- Normalize counts with V4 in vitro counts and log transform
- Very high correlation, DARTFISH BA8 correlates better with DARTFISH PGP1f than RNA-Seq BA8
- Could be the low integer values (many genes with only 1-3 counts) combined with same normalization count in denominator make correlation look better
File:V4 log norm PGP1f vs BA8.jpeg
PGP1f: DARTFISH vs RNA-Seq[edit]
- Data from http://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-54733/
- RNA-Seq data cited by FISSEQ paper
- File: GSM1322928_results_s_6.txt
- Convert RefSeq mRNA ID to HUGO symbol with Biomart
- Some genes have multiple RefSeq IDs for one HUGO symbol
- Use sum of isoform FPKM values for gene FPKM value
- V4
File:PGP1fV4 logDARTFISH vs logRNASeq regression.jpeg
- V7
File:PGP1fV7 logDARTFISH vs logRNASeq regression.jpeg
PGP1f: DARTFISH vs Illumina BeadArray[edit]
- http://genome-tech.ucsd.edu/public/iPSData/Jay_lee_Jan_2008_Illumina_gene_expression.xls
- Used Biomart to convert RefSeq mRNA Acc ID to HGNC symbol
- Many did not have an HGNC symbol but these genes were mostly pseudogenes or other anomalies
- Created file with only HGNC symbol and gene expression (no headers): Jay_lee_Jan_2008_Illumina_gene_expression.txt
- Multiple rows with same gene...
- Take average value: Jay_lee_Jan_2008_Illumina_gene_expression_avgcollapsed.txt
- Take max value: Jay_lee_Jan_2008_Illumina_gene_expression_maxcollapsed.txt
- Multiple rows with same gene...
Using Average[edit]
- V4
File:PGP1fV4 logDARTFISH vs logBeadArrayAvg regression.jpeg
- V7
File:PGP1fV7 logDARTFISH vs logBeadArrayAvg regression.jpeg
Using Max[edit]
- V4
File:PGP1fV4 logDARTFISH vs logBeadArrayMax regression.jpeg
- V7
File:PGP1fV7 logDARTFISH vs logBeadArrayMax regression.jpeg
PGP1f: DARTFISH vs FISSEQv1[edit]
- V4V7_Genes_CompareJayFISSEQData.xlsx
- Counts are from FISSEQ_Science_2014_SeqData_csv\Fibroblast.FISSEQ.2013.06.06.RHex.29bp
- V4
File:PGP1fV4 logDARTFISH vs logFISSEQv1 regression.jpeg
- V7
File:PGP1fV7 logDARTFISH vs logFISSEQv1 regression.jpeg
PGP1f: DARTFISH vs FISSEQv2[edit]
- FISSEQvsCellTypeSpecificRNA_Seq.xlsx
- V4
File:PGP1fV4 logDARTFISH vs logFISSEQv2 regression.jpeg
- V7
File:PGP1fV7 logDARTFISH vs logFISSEQv2 regression.jpeg
PGP1f: DARTFISH vs RNA-Seqv2[edit]
- FISSEQvsCellTypeSpecificRNA_Seq.xlsx
- Has a column labeled "fibro" which I think is RNA-Seq data? Must ask Hosuk where he got it from
- V4
File:PGP1fV4 logDARTFISH vs logRNASeqv2 regression.jpeg
- V7
File:PGP1fV7 logDARTFISH vs logRNASeqv2 regression.jpeg
FISSEQv1 vs FISSEQv2[edit]
- Very few data points because only including genes in V4/V7 probeset AND that are present in both FISSEQv1 and FISSEQv2 data sets