Kun:LabNotes/SingleCellExpr/2011-9-7: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
 
(16 intermediate revisions by the same user not shown)
Line 78: Line 78:
   /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_9_vs_10 tophat_Mm_Indx9/accepted_hits.bam tophat_Mm_Indx10/accepted_hits.bam  
   /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_9_vs_10 tophat_Mm_Indx9/accepted_hits.bam tophat_Mm_Indx10/accepted_hits.bam  
   /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_14_vs_15 tophat_Mm_Indx14/accepted_hits.bam tophat_Mm_Indx15/accepted_hits.bam  
   /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_14_vs_15 tophat_Mm_Indx14/accepted_hits.bam tophat_Mm_Indx15/accepted_hits.bam  
  /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_5_vs_9 tophat_Mm_Indx5/accepted_hits.bam tophat_Mm_Indx9/accepted_hits.bam
  /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_6_vs_11 tophat_Mm_Indx6/accepted_hits.bam tophat_Mm_Indx11/accepted_hits.bam
  /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_11_vs_12 tophat_Mm_Indx11/accepted_hits.bam tophat_Mm_Indx12/accepted_hits.bam
*Then I plotted the FPKM values between each pairs of replicates.
*Then I plotted the FPKM values between each pairs of replicates.
   [[Image:FPKM_Indx3_vs_Indx7.png]][[Image:FPKM_Indx4_vs_Indx8.png]]
  Between two batches of experiments:
   [[Image:FPKM_Indx9_vs_Indx10.png]][[Image:FPKM_Indx14_vs_Indx15.png]]
                R=0.939                                    R=0.928
   [[Image:FPKM_Indx3_vs_Indx7.png|300px]][[Image:FPKM_Indx4_vs_Indx8.png|300px]]
                R=0.935                                    R=0.978
  [[Image:FPKM_Indx5_vs_Indx9.png|300px]][[Image:FPKM_Indx6_vs_Indx11.png|300px]]
 
  Within one batch of experiment.
                R=0.911                                    R=0.984 
   [[Image:FPKM_Indx9_vs_Indx10.png|300px]][[Image:FPKM_Indx11_vs_Indx12.png|300px]]
                R=0.552 (a few transcripts seem to be missing in one sample)
  [[Image:FPKM_Indx14_vs_Indx15.png|300px]]
===Looking for differentially expressed genes===
*I just realized that when feeding cuffdiff with a gtf file generated by cuffcompare, all the gene names are XLOC_xxxxx, which are difficult to link to the known genes. So I changed the gtf file to the UCSC mouse genes.gtf.
 
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_E9.5_vs_E11.5 tophat_Mm_Indx3/accepted_hits.bam,tophat_Mm_Indx7/accepted_hits.bam tophat_Mm_Indx4/accepted_hits.bam,tophat_Mm_Indx8/accepted_hits.bam
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_E11.5_vs_wtE13.5f tophat_Mm_Indx4/accepted_hits.bam,tophat_Mm_Indx8/accepted_hits.bam tophat_Mm_Indx6/accepted_hits.bam,tophat_Mm_Indx11/accepted_hits.bam,tophat_Mm_Indx12/accepted_hits.bam
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_E11.5_vs_wtE13.5m tophat_Mm_Indx4/accepted_hits.bam,tophat_Mm_Indx8/accepted_hits.bam tophat_Mm_Indx5/accepted_hits.bam,tophat_Mm_Indx9/accepted_hits.bam,tophat_Mm_Indx10/accepted_hits.bam
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_wtE13.5m_vs_wtE13.5f tophat_Mm_Indx5/accepted_hits.bam,tophat_Mm_Indx9/accepted_hits.bam,tophat_Mm_Indx10/accepted_hits.bam tophat_Mm_Indx6/accepted_hits.bam,tophat_Mm_Indx11/accepted_hits.bam,tophat_Mm_Indx12/accepted_hits.bam
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_wtE13.5f_vs_koE13.5f tophat_Mm_Indx6/accepted_hits.bam,tophat_Mm_Indx11/accepted_hits.bam,tophat_Mm_Indx12/accepted_hits.bam  tophat_Mm_Indx14/accepted_hits.bam,tophat_Mm_Indx15/accepted_hits.bam
          Wild-type: E9.5 vs. E11.5          Wild-type: E11.5 vs. E13.5-female
  [[Image:FPKM_wtE9.5_vs_wtE11.5.png|300px]][[Image:FPKM_wtE11.5_vs_wtE13.5f.png|300px]]
                  [[Media:FPKM_wtE9.5_vs_wtE11.5.xlsx|Data]]                                      [[Media:FPKM_wtE11.5_vs_wtE13.5f.png.xlsx|Data]]
    Wild-type: E113.5-male vs. E13.5-female    E13.5-female: wild-type vs. TET1-KO
  [[Image:FPKM_wtE13.5m_vs_wtE13.5f.png|300px]][[Image:FPKM_wtE13.5f_vs_koE13.5f.png|300px]]
                  [[Media:FPKM_wtE13.5m_vs_wtE13.5f.xlsx|Data]]                                      [[Media:FPKM_wtE13.5f_vs_koE13.5f.png.xlsx|Data]]
*I randomly picked a few genes that showed biggest differences, and checked their functions. Clearly the data make sense biologically. Examples include:
**Dazl and Sohlh2 were up-regulated by over 100-fold from E9.5 to E11.5.
**Taf7l was up-regulated by ~50 folds from E11.5 to E13.5f.
**Lefty1 and Lefty2 expressed at >100x higher in E13.5m than in E13.5f.
**Tet1's expression was almost completely lost in ko-E13.5f.
===Global analysis of all data sets===
*First I need to combine the gene expression levels (FPKM) for all data sets in one matrix. For this I write a script [[Media:combine_geneTracking.txt|combine_geneTracking.pl]], and ran it in the directory:/home/kunzhang/SingleCell/Data/HL101
*For hierarchical clustering analysis, missing values are allowed in the matrix. I loaded the matrix into R, and create a dendrogram using the following commands:
  x=read.table("Mm_Indx3-15_FPKM_matrix.txt",header=TRUE,row.names=1);
  x.cor = cor(x,use="pairwise.complete.obs",method="pearson")
  x.d = as.dist(1-x.cor)
  x.h = hclust(x.d,"average")
  plot(x.h)
 
  [[Image:PGC-RNAseq-Hierarchical_analysis_092011.png|400px]]
*Replicates tend to cluster together, suggesting that the quality of the data is good.
*For PCA analysis, no missing value is allowed. So I set $no_missing_values = 1, and re-created the matrix. The PCA plot was generated with the following commands:
  methylMatrix=read.table("Mm_Indx3-15_FPKM_matrix_no_NA.txt",header=TRUE,row.names=1)
  pca1 <- prcomp(t(methylMatrix), scale=TRUE)
  mycolors=c("cyan", "green", "blue", "magenta","cyan", "green", "blue", "blue","magenta","magenta","dark  red","red","red");
  plot(pca1$x, pch=20, col=mycolors);
  plot(pca1$x, type="n")
  text(pca1$x, rownames(pca1$x), cex=0.8, col=mycolors)
 
  [[Image:PGC-RNAseq-PCA_plot_092011.png|500px]]
*PC1 seems to explain the difference in developmental stage, whereas PC2 accounts for biological or technical variability among replicates.
*Tet1 knock-out seem to slow down the development into the E13.5 stage.

Latest revision as of 23:58, 26 September 2011

RNAseq data analysis[edit]

Analysis of low-pass sequencing data.[edit]

  • The first round of low-pass sequencing was done in HL101 run (SE60bp+index), lane 7-8.
  • I ran tophat on each data set, using mouse gene annotation as the model.
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx3.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx3.txt.gz > RNAseq_Indx3.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx4.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx4.txt.gz > RNAseq_Indx4.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx5.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx5.txt.gz > RNAseq_Indx5.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx6.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx6.txt.gz > RNAseq_Indx6.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx7.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx7.txt.gz > RNAseq_Indx7.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx8.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx8.txt.gz > RNAseq_Indx8.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx9.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx9.txt.gz > RNAseq_Indx9.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx10.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx10.txt.gz > RNAseq_Indx10.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx11.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx11.txt.gz > RNAseq_Indx11.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx12.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx12.txt.gz > RNAseq_Indx12.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx13.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx13.txt.gz > RNAseq_Indx13.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx14.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx14.txt.gz > RNAseq_Indx14.fq
  less /home/kunzhang/SeqStore/110902_HL101/s_7_Indx15.txt.gz /home/kunzhang/SeqStore/110902_HL101/s_8_Indx15.txt.gz > RNAseq_Indx15.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx3 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx3.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx4 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx4.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx5 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx5.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx6 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx6.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx7 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx7.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx8 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx8.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx9 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx9.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx10 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx10.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx11 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx11.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx12 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx12.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx13 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx13.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx14 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx14.fq
  /home/kunzhang/softwares/tophat-latest/tophat -p 8 --solexa1.3-quals -o tophat_Mm_Indx15 -G /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf  /GenomeDB/UCSC/Mus_musculus/mm9/Sequence/BowtieIndex/genome RNAseq_Indx15.fq
  • Then I checked the tophat_Mm_IndxN/logs/bowtie.left_kept_reads.fixmap.log files, and used the mapping rates in these files as the proxy of the true mapping rates.
Sample Name Index Qualty-filtered reads Uniquely mappable reads Non-mappable reads Mapping rate
E9.5 Indx3 3,775,311 2,715,872 1,014,808 71.9%
E11.5 Indx4 3,987,505 2,710,512 1,238,018 68.0%
wt E13.5m1 Indx5 4,212,162 2,921,411 1,250,730 69.4%
wt E13.5f1 Indx6 3,557,137 2,506,010 1,019,618 70.5%
E9.5 Indx7 5,628,397 3,964,069 1,571,087 70.4%
E11.5 Indx8 6,095,896 4,322,793 1,687,193 70.9%
wt E13.5m1 Indx9 6,284,443 4,290,999 1,925,123 68.3%
wt E13.5m2 Indx10 4,649,460 3,276,710 1,316,390 70.5%
wt E13.5f1 Indx11 5,953,433 4,199,249 1,693,646 70.5%
wt E13.5f2 Indx12 6,712,315 4,671,588 1,969,301 69.6%
KO E13.5m Indx13 4,954,624 3,517,633 1,380,709 71.0%
KO E13.5f1 Indx14 7,223,059 5,137,419 1,999,302 71.1%
KO E13.5f2 Indx15 8,053,410 5,681,254 2,272,034 70.5%
  • Conclusions:
    • The mapping rate is consistently at ~70% across all libraries.
    • To generate roughly 20 millions mappable reads in each library, we need 4-5 HiSeq sequencing lanes.

Further examination of the data quality[edit]

  • I first create a master transcripts.gtf file that combined all transcripts found in the 13 data sets.
 /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffcompare -o cuffcompare_all -R tophat_Mm_Indx3/transcripts.gtf tophat_Mm_Indx4/transcripts.gtf ... tophat_Mm_Indx15/transcripts.gtf 
  • Then I ran cuffdiff on technical/biological replicates to examine the noise in the data.
 /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_3_vs_7 tophat_Mm_Indx3/accepted_hits.bam tophat_Mm_Indx7/accepted_hits.bam 
 /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_4_vs_8 tophat_Mm_Indx4/accepted_hits.bam tophat_Mm_Indx8/accepted_hits.bam 
 /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_9_vs_10 tophat_Mm_Indx9/accepted_hits.bam tophat_Mm_Indx10/accepted_hits.bam 
 /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_14_vs_15 tophat_Mm_Indx14/accepted_hits.bam tophat_Mm_Indx15/accepted_hits.bam 
 /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_5_vs_9 tophat_Mm_Indx5/accepted_hits.bam tophat_Mm_Indx9/accepted_hits.bam 
 /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_6_vs_11 tophat_Mm_Indx6/accepted_hits.bam tophat_Mm_Indx11/accepted_hits.bam 
 /home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff cuffcompare_all.combined.gtf -N -p 8 -o cuffDiff_11_vs_12 tophat_Mm_Indx11/accepted_hits.bam tophat_Mm_Indx12/accepted_hits.bam 
  • Then I plotted the FPKM values between each pairs of replicates.
 Between two batches of experiments:
                R=0.939                                    R=0.928
 File:FPKM Indx3 vs Indx7.pngFile:FPKM Indx4 vs Indx8.png
                R=0.935                                    R=0.978
 File:FPKM Indx5 vs Indx9.pngFile:FPKM Indx6 vs Indx11.png
 Within one batch of experiment.
                R=0.911                                    R=0.984   
 File:FPKM Indx9 vs Indx10.pngFile:FPKM Indx11 vs Indx12.png
                R=0.552 (a few transcripts seem to be missing in one sample)
 File:FPKM Indx14 vs Indx15.png

Looking for differentially expressed genes[edit]

  • I just realized that when feeding cuffdiff with a gtf file generated by cuffcompare, all the gene names are XLOC_xxxxx, which are difficult to link to the known genes. So I changed the gtf file to the UCSC mouse genes.gtf.
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_E9.5_vs_E11.5 tophat_Mm_Indx3/accepted_hits.bam,tophat_Mm_Indx7/accepted_hits.bam tophat_Mm_Indx4/accepted_hits.bam,tophat_Mm_Indx8/accepted_hits.bam 
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_E11.5_vs_wtE13.5f tophat_Mm_Indx4/accepted_hits.bam,tophat_Mm_Indx8/accepted_hits.bam tophat_Mm_Indx6/accepted_hits.bam,tophat_Mm_Indx11/accepted_hits.bam,tophat_Mm_Indx12/accepted_hits.bam
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_E11.5_vs_wtE13.5m tophat_Mm_Indx4/accepted_hits.bam,tophat_Mm_Indx8/accepted_hits.bam tophat_Mm_Indx5/accepted_hits.bam,tophat_Mm_Indx9/accepted_hits.bam,tophat_Mm_Indx10/accepted_hits.bam
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_wtE13.5m_vs_wtE13.5f tophat_Mm_Indx5/accepted_hits.bam,tophat_Mm_Indx9/accepted_hits.bam,tophat_Mm_Indx10/accepted_hits.bam tophat_Mm_Indx6/accepted_hits.bam,tophat_Mm_Indx11/accepted_hits.bam,tophat_Mm_Indx12/accepted_hits.bam
/home/kunzhang/softwares/cufflinks-1.0.3.Linux_x86_64/cuffdiff /GenomeDB/UCSC/Mus_musculus/mm9/Annotation/Genes/genes.gtf -N -p 8 -o cuffDiff_wtE13.5f_vs_koE13.5f tophat_Mm_Indx6/accepted_hits.bam,tophat_Mm_Indx11/accepted_hits.bam,tophat_Mm_Indx12/accepted_hits.bam  tophat_Mm_Indx14/accepted_hits.bam,tophat_Mm_Indx15/accepted_hits.bam
          Wild-type: E9.5 vs. E11.5           Wild-type: E11.5 vs. E13.5-female
  File:FPKM wtE9.5 vs wtE11.5.pngFile:FPKM wtE11.5 vs wtE13.5f.png 
                  Data                                       Data 
   Wild-type: E113.5-male vs. E13.5-female    E13.5-female: wild-type vs. TET1-KO
  File:FPKM wtE13.5m vs wtE13.5f.pngFile:FPKM wtE13.5f vs koE13.5f.png
                  Data                                       Data
  • I randomly picked a few genes that showed biggest differences, and checked their functions. Clearly the data make sense biologically. Examples include:
    • Dazl and Sohlh2 were up-regulated by over 100-fold from E9.5 to E11.5.
    • Taf7l was up-regulated by ~50 folds from E11.5 to E13.5f.
    • Lefty1 and Lefty2 expressed at >100x higher in E13.5m than in E13.5f.
    • Tet1's expression was almost completely lost in ko-E13.5f.

Global analysis of all data sets[edit]

  • First I need to combine the gene expression levels (FPKM) for all data sets in one matrix. For this I write a script combine_geneTracking.pl, and ran it in the directory:/home/kunzhang/SingleCell/Data/HL101
  • For hierarchical clustering analysis, missing values are allowed in the matrix. I loaded the matrix into R, and create a dendrogram using the following commands:
  x=read.table("Mm_Indx3-15_FPKM_matrix.txt",header=TRUE,row.names=1);
  x.cor = cor(x,use="pairwise.complete.obs",method="pearson")
  x.d = as.dist(1-x.cor)
  x.h = hclust(x.d,"average")
  plot(x.h)
  File:PGC-RNAseq-Hierarchical analysis 092011.png
  • Replicates tend to cluster together, suggesting that the quality of the data is good.
  • For PCA analysis, no missing value is allowed. So I set $no_missing_values = 1, and re-created the matrix. The PCA plot was generated with the following commands:
  methylMatrix=read.table("Mm_Indx3-15_FPKM_matrix_no_NA.txt",header=TRUE,row.names=1)
  pca1 <- prcomp(t(methylMatrix), scale=TRUE)
  mycolors=c("cyan", "green", "blue", "magenta","cyan", "green", "blue", "blue","magenta","magenta","dark   red","red","red");
  plot(pca1$x, pch=20, col=mycolors); 
  plot(pca1$x, type="n") 
  text(pca1$x, rownames(pca1$x), cex=0.8, col=mycolors)
  File:PGC-RNAseq-PCA plot 092011.png
  • PC1 seems to explain the difference in developmental stage, whereas PC2 accounts for biological or technical variability among replicates.
  • Tet1 knock-out seem to slow down the development into the E13.5 stage.