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
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 #
|
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
'
|
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
'
|
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
File:CD 496.jpg
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