Kun:LabNotes/SingleCellExpr/2014-6-5
Jump to navigation
Jump to search
Preparing the 2014Jun data set
Data source
- The data came from 18 C1 experiments.
- 20131219_C1_hNuc-20140222_Expt123Hiseq(oligo dT only)
- 20140123_C1_hNuc-20140210_Expt119HiSeq
- 20140205_C1_hNuc-20140222_Expt123Hiseq
- 20140212_C1_hNuc-20140222_Expt123Hiseq
- 20140213_C1_hNuc-20140222_Expt123Hiseq
- 20140218_C1_hNuc-20140317_Expt128HiSeq
- 20140220_C1_hNuc-20140317_Expt128HiSeq
- 20140225_C1_hNuc-20140317_Expt128HiSeq
- 20140226_C1_hNuc-20140317_Expt128HiSeq
- 20140227_C1_hNuc-20140414_Expt131HiSeq
- 20140313_C1_hNuc-20140414_Expt131HiSeq
- 20140317_C1_hNuc-20140414_Expt131HiSeq
- 20140304_C1_hNuc-20140521_Expt138HiSeq
- 20140318_C1_hNuc-20140521_Expt138HiSeq
- 20140320_C1_hNuc-20140521_Expt138HiSeq
- 20140321_C1_hNuc-20140521_Expt138HiSeq
cd /home/kunzhang/RNAseq/SCAP/expr_matrix ../scripts/merge_expr_matrix.pl 20131219_C1_hNuc-20140222_TPM_matrix_name_corrected.txt 20140210_96Samples_hNuc_Expt119HiSeq_C1_STAR_0123_TPM_matrix_name_corrected.txt 20140205_C1_hNuc-20140222_TPM_matrix_name_corrected.txt 20140212_C1_hNuc-20140222_TPM_matrix_name_corrected.txt 20140213_C1_hNuc-20140222_TPM_matrix_name_corrected.txt 20140218_C1_hNuc_expr_matrix.txt 20140220_C1_hNuc_expr_matrix.txt 20140225_C1_hNuc_expr_matrix.txt 20140226_C1_hNuc_expr_matrix.txt 20140227_C1_hNuc_expr_matrix.txt 20140311_C1_hNuc_expr_matrix.txt 20140313_C1_hNuc_expr_matrix.txt 20140317_C1_hNuc_expr_matrix.txt 20140304_C1_hNuc_expr_matrix.txt 20140318_C1_hNuc_expr_matrix.txt 20140320_C1_hNuc_expr_matrix.txt 20140321_C1_hNuc_expr_matrix.txt > 2014_1st_Quater_18_expr_matrix.txt cd ../summary_stats ../scripts/get_excluded_sample_id.pl 20131219_C1_hNuc-20140222_mappnig_stats_name_corrected.txt 20140123_C1_hNuc-20140210_mappnig_stats_name_corrected.txt 20140205_C1_hNuc-20140222_mapping_stats_name_corrected.txt 20140212_C1_hNuc-20140222_mapping_stats_name_corrected.txt 20140212_C1_hNuc-20140222_mapping_stats.txt 20140213_C1_hNuc-20140222_mapping_stats_name_corrected.txt 20140218_C1_mapping_stat.txt 20140220_C1_mapping_stat.txt 20140225_C1_mapping_stat.txt 20140226_C1_mapping_stat.txt 20140227_C1_hNuc_mapping_stats.txt 20140311_C1_hNuc_mapping_stats.txt 20140313_C1_hNuc_mapping_stats.txt 20140317_C1_hNuc_mapping_stats.txt 20140304_C1_mapping_stat.txt 20140318_C1_mapping_stat.txt 20140320_C1_mapping_stat.txt 20140321_C1_mapping_stat.txt > 2014_1st_Quater_18_expr_matrix_1N_excluded_ids.txt ../scripts/filterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt ../summary_stats/2014_1st_Quater_18_expr_matrix_1N_excluded_ids.txt > 2014_1st_Quater_18_1N_expr_matrix_filtered.txt
A total of 728 data sets (10323 genes) were obtained after the filtering (at least 1M total reads; mapping rate >0.4; # ERCC detected > 24; % ERCC reads < 0.3; <30% missing values per gene).
Iterative over-dispersion analysis
Round 1
Identify over-dispersed genes
#Find over-dispersed genes for the entire data set, using a Z cutoff of 2. ../scripts/exprMatrix2overDispersedGenes.pl 2014_1st_Quater_18_1N_expr_matrix_filtered.txt 2 > 2014_1st_Quater_18_1N_Z2_overDispersed_genes.txt #Extract the gene expression matrix for the over-dispersed genes. ../scripts/extract_gene_set_ENSG.pl 2014_1st_Quater_18_1N_Z2_overDispersed_genes.txt 2014_1st_Quater_18_1N_expr_matrix_filtered.txt > 2014_1st_Quater_18_1N_Z2_overDispersed_expr_matrix.txt
Hierarchical clustering of samples based on these genes
#Visualize the sample clustering. library(gplots) x=read.table("2014_1st_Quater_18_1N_Z2_overDispersed_expr_matrix.txt",header=TRUE,row.names=1); x.cor = cor(log10(x+1),use="pairwise.complete.obs",method="pearson") heatmap.2(as.matrix(x.cor), col=bluered(128), scale="none", cexCol=0.1, cexRow=0.1 ,key=T, symkey=F,density.info="histogram",trace="none",dendrogram="both",Rowv=TRUE,Colv=TRUE) File:2014 1st Quater 18 over-dispersed-Z2 sample clustering.png #Cut the dendrogram, and extract the sample names in each cluster hc.rows<- hclust(dist(x.cor)) plot(hc.rows,cex=0.1) ct<- cutree(hc.rows, h=11) # the threshold (11) was selected by visually inspect the dendrogram. rect.hclust(hc.rows, h=11) # draw red rectangles to mark the subgroups write.table(ct, file="2014_1st_Quater_18_1N_Z2_overDispersed_cutree_3_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE) File:2014 1st Quater 18 over-dispersed-Z2 dendrogram cutree h11.png #Generate gene expression matrix for each cluster ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z2_overDispersed_cutree_3_clusters.txt 1 > 2014_1st_Quater_18_expr_matrix_cutree-h11_c1.txt ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z2_overDispersed_cutree_3_clusters.txt 2 > 2014_1st_Quater_18_expr_matrix_cutree-h11_c2.txt & ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z2_overDispersed_cutree_3_clusters.txt 3 > 2014_1st_Quater_18_expr_matrix_cutree-h11_c3.txt & #Check whether samples in these cluster have specific neural transmitters. ./make_gene_set_heatmap_Hs.pl neural_transmitter_metabolism.txt ../expr_matrix/2014_1st_Quater_18_expr_matrix_cutree-h11_c1.txt ./make_gene_set_heatmap_Hs.pl neural_transmitter_metabolism.txt ../expr_matrix/2014_1st_Quater_18_expr_matrix_cutree-h11_c2.txt ./make_gene_set_heatmap_Hs.pl neural_transmitter_metabolism.txt ../expr_matrix/2014_1st_Quater_18_expr_matrix_cutree-h11_c3.txt Cluster #1 Cluster #2 Cluster #3 File:2014 1st Quater 18 expr matrix cutree-h11 c1 neural transmitter metabolism expr heatmap.pngFile:2014 1st Quater 18 expr matrix cutree-h11 c2 neural transmitter metabolism expr heatmap.pngFile:2014 1st Quater 18 expr matrix cutree-h11 c3 neural transmitter metabolism expr heatmap.png
- Clear Cluster #1 are mostly GAD1-/GAD2- neurons, and Cluster #3 are mostly GAD1+/GAD2+ neurons. Cluster #2 is a mixture of both. While overall there is a separate between excitatory versus inhibitory neurons, the clustering is not exactly based on such a difference. There is probably another feature that determines the clustering.
K-mean partitioning based on the over-dispersed genes
library(cluster) library(fpc) x=read.table("2014_1st_Quater_18_1N_Z2_overDispersed_expr_matrix.txt",header=TRUE,row.names=1) x[is.na(x)] <- 0 y=log10(x+1) tys = scale(t(y)) k_5=kmeans(tys, 5) k_4=kmeans(tys, 4) k_3=kmeans(tys, 3) k_2=kmeans(tys, 2) plotcluster(tys,k_2$cluster) plotcluster(tys,k_3$cluster) plotcluster(tys,k_4$cluster) plotcluster(tys,k_5$cluster) write.table(k_3$cluster, file="2014_1st_Quater_18_1N_Z2_overDispersed_k-mean_3_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE) write.table(k_2$cluster, file="2014_1st_Quater_18_1N_Z2_overDispersed_k-mean_2_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE) File:2014 1st Quater 18 over-dispersed-Z2 k2 clusters.pngFile:2014 1st Quater 18 over-dispersed-Z2 k3 clusters.png File:2014 1st Quater 18 over-dispersed-Z2 k4 clusters.pngFile:2014 1st Quater 18 over-dispersed-Z2 k5 clusters.png #Generate gene expression matrix for each cluster ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z2_overDispersed_k-mean_2_clusters.txt 1 > 2014_1st_Quater_18_expr_matrix_k-mean2_c1.txt & ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z2_overDispersed_k-mean_2_clusters.txt 2 > 2014_1st_Quater_18_expr_matrix_k-mean2_c2.txt & ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z2_overDispersed_k-mean_3_clusters.txt 1 > 2014_1st_Quater_18_expr_matrix_k-mean3_c1.txt & ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z2_overDispersed_k-mean_3_clusters.txt 2 > 2014_1st_Quater_18_expr_matrix_k-mean3_c2.txt & ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z2_overDispersed_k-mean_3_clusters.txt 3 > 2014_1st_Quater_18_expr_matrix_k-mean3_c3.txt &
./make_gene_set_heatmap_Hs.pl neural_transmitter_metabolism.txt ../expr_matrix/2014_1st_Quater_18_expr_matrix_k-mean2_c1.txt ./make_gene_set_heatmap_Hs.pl neural_transmitter_metabolism.txt ../expr_matrix/2014_1st_Quater_18_expr_matrix_k-mean2_c2.txt ./make_gene_set_heatmap_Hs.pl neural_transmitter_metabolism.txt ../expr_matrix/2014_1st_Quater_18_expr_matrix_k-mean3_c1.txt ./make_gene_set_heatmap_Hs.pl neural_transmitter_metabolism.txt ../expr_matrix/2014_1st_Quater_18_expr_matrix_k-mean3_c2.txt ./make_gene_set_heatmap_Hs.pl neural_transmitter_metabolism.txt ../expr_matrix/2014_1st_Quater_18_expr_matrix_k-mean3_c3.txt Cluster #1 (k=2) Cluster #2 (k=2) File:2014 1st Quater 18 expr matrix k-mean2 c1 neural transmitter metabolism expr heatmap.pngFile:2014 1st Quater 18 expr matrix k-mean2 c2 neural transmitter metabolism expr heatmap.png Cluster #1(k=3) Cluster #2(k=3) Cluster #3(k=3) File:2014 1st Quater 18 expr matrix k-mean3 c1 neural transmitter metabolism expr heatmap.pngFile:2014 1st Quater 18 expr matrix k-mean3 c2 neural transmitter metabolism expr heatmap.pngFile:2014 1st Quater 18 expr matrix k-mean3 c3 neural transmitter metabolism expr heatmap.png
- The sample clusters obtained by cutting the dendrogram are not entirely consistent with the k-mean generated clusters. For example, Cluster #1 in the k-3 clustering result include all #3, plus some #2 samples in the dendrogram. Similarly, GAD1/GAD2 seem to contribute to the separation, but there is also something else.
Identify over-dispersed genes using a less stringent threshold
- When looked into the matrix of over-dispersed genes, I found that only 95 genes were in the list. Perhaps that's too few and not sufficient for robust sample clustering. Therefore I repeated the whole analysis using a Z cutoff of 1.
../scripts/exprMatrix2overDispersedGenes.pl 2014_1st_Quater_18_1N_expr_matrix_filtered.txt 1 > 2014_1st_Quater_18_1N_Z1_overDispersed_genes.txt #Extract the gene expression matrix for the over-dispersed genes. ../scripts/extract_gene_set_ENSG.pl 2014_1st_Quater_18_1N_Z1_overDispersed_genes.txt 2014_1st_Quater_18_1N_expr_matrix_filtered.txt > 2014_1st_Quater_18_1N_Z1_overDispersed_expr_matrix.txt
x=read.table("2014_1st_Quater_18_1N_Z1_overDispersed_expr_matrix.txt",header=TRUE,row.names=1); x.cor = cor(log10(x+1),use="pairwise.complete.obs",method="pearson") hm<-heatmap.2(as.matrix(x.cor), col=bluered(128), scale="none", cexCol=0.1, cexRow=0.1 ,key=T, symkey=F,density.info="histogram",trace="none",dendrogram="both",Rowv=TRUE,Colv=TRUE)
hc.rows<- hclust(dist(x.cor)) plot(hc.rows,cex=0.1) ct<- cutree(hc.rows, h=11) # rect.hclust(hc.rows, h=11) # draw red rectangles to mark the subgroups write.table(ct, file="2014_1st_Quater_18_1N_Z1_overDispersed_cutree_2_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE)
x=read.table("2014_1st_Quater_18_1N_Z1_overDispersed_expr_matrix.txt",header=TRUE,row.names=1) x[is.na(x)] <- 0 y=log10(x+1) tys = scale(t(y)) k_5=kmeans(tys, 5) k_4=kmeans(tys, 4) k_3=kmeans(tys, 3) k_2=kmeans(tys, 2) plotcluster(tys,k_2$cluster) plotcluster(tys,k_3$cluster) plotcluster(tys,k_4$cluster) plotcluster(tys,k_5$cluster) write.table(k_3$cluster, file="2014_1st_Quater_18_1N_Z1_overDispersed_k-mean_3_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE) write.table(k_2$cluster, file="2014_1st_Quater_18_1N_Z1_overDispersed_k-mean_2_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE) ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z1_overDispersed_cutree_2_clusters.txt 1 > 2014_1st_Quater_18_expr_matrix_Z1_cutree-h11_c1.txt ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z1_overDispersed_cutree_2_clusters.txt 2 > 2014_1st_Quater_18_expr_matrix_Z1_cutree-h11_c2.txt ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z1_overDispersed_k-mean_2_clusters.txt 1 > 2014_1st_Quater_18_expr_matrix_Z1_k-mean2_c1.txt ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z1_overDispersed_k-mean_2_clusters.txt 2 > 2014_1st_Quater_18_expr_matrix_Z1_k-mean2_c2.txt ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z1_overDispersed_k-mean_3_clusters.txt 1 > 2014_1st_Quater_18_expr_matrix_Z1_k-mean3_c1.txt ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z1_overDispersed_k-mean_3_clusters.txt 2 > 2014_1st_Quater_18_expr_matrix_Z1_k-mean3_c2.txt ../scripts/extractSampleClusterExprMatrix.pl 2014_1st_Quater_18_expr_matrix.txt 2014_1st_Quater_18_1N_Z1_overDispersed_k-mean_3_clusters.txt 3 > 2014_1st_Quater_18_expr_matrix_Z1_k-mean3_c3.txt