Kun:LabNotes/SingleCellExpr/2011-9-7: Difference between revisions
Jump to navigation
Jump to search
(14 intermediate revisions by the same user not shown) | |||
Line 82: | Line 82: | ||
/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 | /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| | Between two batches of experiments: | ||
[[Image:FPKM_Indx9_vs_Indx10.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]
- The data came from limited number of mouse PGC cells provided by Yi Zhang's lab. Rui made the RNA sequencing libraries.
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.