Rui:LabNotes/SingleCell/2014-7-4

From ZhangLabWiki
Jump to navigation Jump to search

Cell type analysis on 20 chips[edit]

  • Goal: general workflow -> cell groups -> signature genes -> FISH candidate
  • Refined analysis in Wei group: noise cancellation + clustering analysis
  • Candidate validation in Jerold/Kun group: FISH (Yun/Blue/Daniel)
  • Additional: bulk sample vs singles

20 chips primary analysis[edit]

  • Primary analysis [1]
  • Data: genome-miner [/media/LTS_33T/RL_LTS33T]
  • Analysis: genome-miner [/home/ruiliu/Ext12T/7.3.14_20chips_hNuc]

SuperMatrix[edit]

  • Quick & dirty assembly of TPM matrix
  • Keys: rename and ensure IDs/numbers consistent in all files
  • Need the noise normalization based on ERCC from Wei's group to refine TPM calling
  • Need to implement more systematic quality filters (literature)
/home/kunzhang/RNAseq/SCAP/scripts/merge_expr_matrix.pl /media/LTS_33T/RL_LTS33T/20131219_C1_hNuc-20140222_Expt123Hiseq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140205_C1_hNuc-20140222_Expt123Hiseq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140123_C1_hNuc-20140210_Expt119HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140212_C1_hNuc-20140222_Expt123Hiseq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140213_C1_hNuc-20140222_Expt123Hiseq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140218_C1_hNuc-20140317_Expt128HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140220_C1_hNuc-20140317_Expt128HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140225_C1_hNuc-20140317_Expt128HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140226_C1_hNuc-20140317_Expt128HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140227_C1_hNuc-20140414_Expt131HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140304_C1_hNuc-20140521_Expt138HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140311_C1_hNuc-20140414_Expt131HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140313_C1_hNuc-20140414_Expt131HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140317_C1_hNuc-20140414_Expt131HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140318_C1_hNuc-20140521_Expt138HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140320_C1_hNuc-20140521_Expt138HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140321_C1_hNuc-20140521_Expt138HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140326_C1_hNuc-20140620_Expt141HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140331_C1_hNuc-20140620_Expt141HiSeq/TPMmatrix.txt /media/LTS_33T/RL_LTS33T/20140402_C1_hNuc-20140620_Expt141HiSeq/TPMmatrix.txt > 20chip.TPMmatrix.txt
./get_matrix_size.pl < 20chip.TPMmatrix.newID-2.txt
  • Simple filter by at least 1M total reads; mapping rate >0.4 => relative_hg19map_rate>0.4/0.2; # ERCC detected > 24 => 22 in fpkm; % ERCC reads < 0.3 => off
./get_excluded_sample_id_hg19-0.4.pl 20chip.mapStat.newID.txt > 20chip.TPMmatrix_hg19-0.4_1N_excluded_ids.txt
./get_excluded_sample_id_hg19-0.2.pl 20chip.mapStat.newID.txt > 20chip.TPMmatrix_hg19-0.2_1N_excluded_ids.txt
./filterExprMatrix-70.pl 20chip.TPMmatrix.newID-2.txt  20chip.TPMmatrix_hg19-0.4_1N_excluded_ids.txt > 20chip.TPMmatrix.newID_0.4filtered_0.7allowed.txt
./filterExprMatrix-70.pl 20chip.TPMmatrix.newID-2.txt  20chip.TPMmatrix_hg19-0.2_1N_excluded_ids.txt > 20chip.TPMmatrix.newID_0.2filtered_0.7allowed.txt
  • Some numbers:
map% singles missing genes filtered TPMmatrix file '
map>0.4 728 0.3 10323
r.hg19>0.4 911 0.5 13333 20chip.TPMmatrix.newID_0.4filtered_0.5allowed.txt
0.7 16978 20chip.TPMmatrix.newID_0.4filtered_0.7allowed.txt
r.hg19>0.2 986 0.5 13148 20chip.TPMmatrix.newID_0.2filtered_0.5allowed.txt
0.7 16831 20chip.TPMmatrix.newID_0.2filtered_0.7allowed.txt

986 singles (16831 genes) passing the criteria of 20% relative hg19 percentage and 70% gene allowed missing[edit]

  • Figures reorganized on [2]

First layer of sample clusters based on overdispersed genes[edit]

Overdispersed genes[edit]

  • 20chip.TPMmatrix.newID_0.2filtered_0.7allowed.txt => 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_genes.txt => 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_TPMmatrix.txt
/home/kunzhang/RNAseq/SCAP/scripts/exprMatrix2overDispersedGenes.pl 20chip.TPMmatrix.newID_0.2filtered_0.7allowed.txt 1 > 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_genes.txt
/home/kunzhang/RNAseq/SCAP/scripts/extract_gene_set_ENSG.pl  20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_genes.txt 20chip.TPMmatrix.newID_0.2filtered_0.7allowed.txt > 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_TPMmatrix.txt
./get_matrix_size.pl < 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_TPMmatrix.txt
#data column: 986
# data line: 557

Hierarchical clustering[edit]

  • hierarchical clustering based on 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_TPMmatrix.txt
 library(gplots)
 x=read.table("20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_TPMmatrix.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:HC 557.jpg
  • However, when I tried simply cluster dendrogram, the cluster makes more sense to me.
x.d = as.dist(1-x.cor)	
x.h = hclust(x.d,"complete")	
plot(x.h, cex = 0.1, cex.main = 2)
File:CD 557.jpg
  • Refine the heatmap2 code by directly cluster dendrogram value without doing further clustering
hc.rows<- hclust(as.dist(1 - x.cor), method="complete")
colv <- as.dendrogram(hc.rows)
rowv <- colv
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=rowv,Colv=colv)
  • Cut dendrogram tree branches
plot(hc.rows,cex=0.1) 									
ct<- cutree(hc.rows, h=0.9) # 									
rect.hclust(hc.rows, h=0.9) # draw red rectangles to mark the subgroups									
write.table(ct, file="xxxx_cutree_2_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE)									
File:HConCD 557.jpg

K means clustering[edit]

  • K-mean clustering based on 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_TPMmatrix.txt
library(cluster) 
library(fpc)
x=read.table("20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_TPMmatrix.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="20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_k-mean_3_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE)
write.table(k_2$cluster, file="20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_k-mean_2_clusters.txt", sep = "\t", row.names=TRUE, col.names=TRUE)
File:557 2c.jpg File:557 3c.jpg 
File:557 4c.jpg File:557 5c.jpg

Comparison b/t two clustering[edit]

Total: 986
2 clusters: 805 samples are same in both method; 181 samples different
3 clusters: 563 are same; 423 different
' 1 2 3 Total
cutree2 737 249 0 986
cutree3 737 144 105 986
kmean2 648 338 0 986
kmean3 563 245 178 986

DEx genes based on sample clustering[edit]

  • By looking at the numbers/plots, hierarchical clustering seems right to me. But as >97% samples are overlapped in both method if only cut into 2 groups, I just stick to cut2 and look for the overlapped DEx genes.
  • I also focus on hierarchical clustering genes too.
  • DEx gene is identified by DESeq2
 library('DESeq2')
 x=read.table("20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_TPMmatrix.txt",header=TRUE,row.names=1);
 colData=read.table("20chip.TPMmatrix_0.2-0.7-Z1_overDispersed_clusters.txt",header=TRUE,row.names=1);
 x[is.na(x)] <- 0
 dds <- DESeqDataSetFromMatrix(countData = round(x), colData = colData, design =~ kmean2)
 dds <- DESeq(dds)
 res <- results(dds)
 resOrdered <- res[order(res$padj),] 
 head(resOrdered,20)
 resSig <- res[ which(res$padj<0.01),]
 head(resSig[order(resSig$log2FoldChange ),], 10)
 tail(resSig[order(resSig$log2FoldChange ),], 10)
  • Gene numbers
' gene #
cutree2 406
kmean2 371
common 339


File:20chips 0.2-0.7-Z1 cutree2.png File:20chips 0.2-0.7-Z1 cut2-both.png

Media:20chips_0.2-0.7-z1_both_DE-0.01.xlsx

Blue's in-depth research on these top 27 genes [3]

top27	log2 fold change
NEAT1	2.003211939
GRIN3A	2.033951837
VAV2	2.062726013
PTPRM	2.080714217
FGF13	2.094894592
RBMS3	2.106795779
MYO16	2.151562063
RP11-123O10.4	2.169064516
KCNIP1	2.200186983
ADRA1A	2.28142167
FAT1	2.321830687
GRIK1	2.363408707
ERBB4	2.368557498
BTBD11	2.467137982
ARX	2.474737422
CTC-321K16.1	2.627130756
MTSS1	2.63965966
BACH1	2.677262451
CXCL14	2.793804062
DNER	2.873455079
ADARB2	2.890106589
GAD1	2.974213273
DLX6-AS1	2.977317694
GAD2	3.233591271
SLC6A1	3.383217471
MAF	3.415295366
SLC6A1-AS1	3.890308415

Second layer of sample clusters based on overdispersed genes[edit]

  • The idea is based on series of bi/tri-clusters
  • The 2 clusters separated from either hierarchical/k-means clustering at the first round are further used to extract sample TPMmatrix followed by overdispersed gene identification/extraction.
  • The 2nd filtered overDispersed geneTPM is used to repeat clustering again

Step1. To generate separate TPM matrix for each group (cluster 1 & 2) based on the first round clustering[edit]

  • . change pct_missing_allowed =70 (70% allowed) from original pl (/home/kunzhang/RNAseq/SCAP/scripts /extractSampleClusterExprMatrix.pl)
./extractSampleClusterExprMatrix_70.pl					
  • . extract sampls of each cluster from original matrix
./extractSampleClusterExprMatrix_70.pl 20chip.TPMmatrix.newID-2.txt 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_cutree_2.txt 1 > 20chip.TPMmatrix_0.2-0.7-Z1_cutree2_c1.txt					
./extractSampleClusterExprMatrix_70.pl 20chip.TPMmatrix.newID-2.txt 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_cutree_2.txt 2 > 20chip.TPMmatrix_0.2-0.7-Z1_cutree2_c2.txt					
./extractSampleClusterExprMatrix_70.pl 20chip.TPMmatrix.newID-2.txt 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_kmean_2.txt 1 > 20chip.TPMmatrix_0.2-0.7-Z1_kmean2_c1.txt					
./extractSampleClusterExprMatrix_70.pl 20chip.TPMmatrix.newID-2.txt 20chip.TPMmatrix.newID_0.2filtered_0.7allowed_Z1_overDispersed_kmean_2.txt 2 > 20chip.TPMmatrix_0.2-0.7-Z1_kmean2_c2.txt					
  • . output table
' clusters samples genes filtered TPMmatrix file
cutree c1 737 17071 20chip.TPMmatrix_0.2-0.7-Z1_cutree2_c1.txt
c2 249 16029 20chip.TPMmatrix_0.2-0.7-Z1_cutree2_c2.txt
kmean c1 648 17482 20chip.TPMmatrix_0.2-0.7-Z1_kmean2_c1.txt
c2 338 15265 20chip.TPMmatrix_0.2-0.7-Z1_kmean2_c2.txt

Step 2: To identify overdispersed genes from whole gene list[edit]

  • . To identify overdispersed genes: [Usage] exprMatrix2overDispersedGenes.pl (filtered_matrix.txt) (z_score_cutoff) > (overdispersed_genelist.txt)
/home/kunzhang/RNAseq/SCAP/scripts/exprMatrix2overDispersedGenes.pl 20chip.TPMmatrix_0.2-0.7-Z1_cutree2_c1.txt 1 > 20chip.TPMmatrix_0.2-0.7-Z1_cutree2_c1_Z1_overDispersed_genes.txt
  • . To refine Ex matrix of overdispersed genes: [Usage] extract_gene_set_ENSG.pl (overdispersed_genelist.txt) (filtered_matrix.txt) > (overdispersed_matrix.txt)
/home/kunzhang/RNAseq/SCAP/scripts/extract_gene_set_ENSG.pl  20chip.TPMmatrix_0.2-0.7-Z1_cutree2_c1_Z1_overDispersed_genes.txt 20chip.TPMmatrix_0.2-0.7-Z1_cutree2_c1.txt > 20chip.TPMmatrix_0.2-0.7-Z1_cutree2-c1-Z1_overDispersed_TPMmatrix.txt	
  • . Output matrix
' clusters samples genes overDisp. '
cutree c1 737 17071 547 20chip.TPMmatrix_0.2-0.7-Z1_cutree2-c1-Z1_overDispersed_TPMmatrix.txt
c2 249 16029 452 20chip.TPMmatrix_0.2-0.7-Z1_cutree2-c2-Z1_overDispersed_TPMmatrix.txt
kmean c1 648 17482 567 20chip.TPMmatrix_0.2-0.7-Z1_kmean2-c1-Z1_overDispersed_TPMmatrix.txt
c2 338 15265 415 20chip.TPMmatrix_0.2-0.7-Z1_kmean2-c2-Z1_overDispersed_TPMmatrix.txt

Step 3: Sample clustering by hierarchical/k-means clustering[edit]

20chip.TPMmatrix_0.2-0.7-Z1_cutree2-c1-Z1_overDispersed_TPMmatrix.txt[edit]
  • This cluster is likely to be excitatory neurons (GAD1-/GAD2-/SLC6A1-)
  • 1st layer cut by cutree2, 2nd layer cut by either cutree or kmeans

File:20chip.TPMmatrix 0.2-0.7-Z1 cutree2-c1-Z1 overDispersed TPMmatrix cutree.jpg

File:20chip.TPMmatrix 0.2-0.7-Z1 cutree2-c1-Z1 overDispersed TPMmatrix kmean2-5.jpg


20chip.TPMmatrix_0.2-0.7-Z1_kmean2-c1-Z1_overDispersed_TPMmatrix.txt[edit]
  • This cluster is likely to be excitatory neurons (GAD1-/GAD2-/SLC6A1-)
  • 1st layer cut by kmean2, 2nd layer cut by either cutree or kmeans

File:20chip.TPMmatrix 0.2-0.7-Z1 kmean2-c1-Z1 overDispersed TPMmatrix cutree.jpg

File:20chip.TPMmatrix 0.2-0.7-Z1 kmean2-c1-Z1 overDispersed TPMmatrix kmean2-5.jpg

911 singles (13333 genes) passing the criteria of 40% relative hg19 percentage and 50% gene allowed missing[edit]

Overdispersed genes[edit]

  • 20chip.TPMmatrix.newID_0.4filtered_0.5allowed.txt => 20chip.TPMmatrix.newID_0.4filtered_0.5allowed_Z1_overDispersed_genes.txt => 20chip.TPMmatrix.newID_0.4filtered_0.5allowed_Z1_overDispersed_TPMmatrix.txt
#data column: 911
# data line: 496

Hierarchical clustering[edit]

  • hierarchical clustering based on 20chip.TPMmatrix.newID_0.4filtered_0.5allowed_Z1_overDispersed_TPMmatrix.txt
File:HC 496.jpg
  • CD figure
File:CD 496.jpg
  • refined heatmap2
File:HC-CD 496.jpg

K means clustering[edit]

  • K-mean clustering based on 20chip.TPMmatrix.newID_0.4filtered_0.5allowed_Z1_overDispersed_TPMmatrix.txt
File:496 2c.jpg File:496 3c.jpg 
File:496 4c.jpg File:496 5c.jpg