Matt:LabNotes/2015-4-8

From ZhangLabWiki
Jump to navigation Jump to search

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]

CountReadsPer_Gene_Probe.pl

CountReadsPer_Gene_Probe_excludeSoftMask.pl

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

File:V4 log PGP1f vs BA8.jpeg

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

File:PGP1f FISSEQv1 vs FISSEQv2 regression.jpeg