Editing
Rui:LabNotes/SingleCell/2014-7-4
Jump to navigation
Jump to search
Warning:
You are not logged in. Your IP address will be publicly visible if you make any edits. If you
log in
or
create an account
, your edits will be attributed to your username, along with other benefits.
Anti-spam check. Do
not
fill this in!
==Cell type analysis on 20 chips== * 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=== * Primary analysis [http://genome-tech.ucsd.edu/LabNotes/index.php/Rui:LabNotes/SingleCell/2014-7-3] * Data: genome-miner [/media/LTS_33T/RL_LTS33T] * Analysis: genome-miner [/home/ruiliu/Ext12T/7.3.14_20chips_hNuc] ===SuperMatrix=== * 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: {| {{table}} class = wikitable | align="center" style="background:#f0f0f0;"|'''map%''' | align="center" style="background:#f0f0f0;"|'''singles''' | align="center" style="background:#f0f0f0;"|'''missing''' | align="center" style="background:#f0f0f0;"|'''genes''' | align="center" style="background:#f0f0f0;"|'''filtered TPMmatrix file''' | align="center" style="background:#f0f0f0;"|'''''' |- | 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== * Figures reorganized on [http://genome-tech.ucsd.edu/LabNotes/index.php/Rui:LabNotes/SingleCell/2014-7-14] ===First layer of sample clusters based on overdispersed genes=== ====Overdispersed genes==== * 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==== * 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) [[Image:HC_557.jpg|800px]] * 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) [[Image:CD_557.jpg|800px]] * 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) [[Image:HConCD_557.jpg|800px]] ====K means clustering==== * 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|400px]] [[File:557_3c.jpg|400px]] [[File:557_4c.jpg|400px]] [[File:557_5c.jpg|400px]] ====Comparison b/t two clustering==== Total: 986 2 clusters: 805 samples are same in both method; 181 samples different 3 clusters: 563 are same; 423 different {| {{table}} | align="center" style="background:#f0f0f0;"|'''''' | align="center" style="background:#f0f0f0;"|'''1''' | align="center" style="background:#f0f0f0;"|'''2''' | align="center" style="background:#f0f0f0;"|'''3''' | align="center" style="background:#f0f0f0;"|'''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==== * 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 {| {{table}} | align="center" style="background:#f0f0f0;"|'''''' | align="center" style="background:#f0f0f0;"|'''gene #''' |- | cutree2||406 |- | kmean2||371 |- | common||339 |} [[File:20chips_0.2-0.7-Z1_cutree2.png|400px]] [[File:20chips_0.2-0.7-Z1_cut2-both.png|400px]] [[Media:20chips_0.2-0.7-z1_both_DE-0.01.xlsx]] Blue's in-depth research on these top 27 genes [http://genome-tech.ucsd.edu/LabNotes/index.php/Blue:RNA-Seq_Analyses:C1_hNuc_1000_Dataset] 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=== * 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==== *. 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 {| {{table}} | align="center" style="background:#f0f0f0;"|'''''' | align="center" style="background:#f0f0f0;"|'''clusters''' | align="center" style="background:#f0f0f0;"|'''samples''' | align="center" style="background:#f0f0f0;"|'''genes''' | align="center" style="background:#f0f0f0;"|'''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==== *. 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 {| {{table}} | align="center" style="background:#f0f0f0;"|'''''' | align="center" style="background:#f0f0f0;"|'''clusters''' | align="center" style="background:#f0f0f0;"|'''samples''' | align="center" style="background:#f0f0f0;"|'''genes''' | align="center" style="background:#f0f0f0;"|'''overDisp.''' | align="center" style="background:#f0f0f0;"|'''''' |- | 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==== =====20chip.TPMmatrix_0.2-0.7-Z1_cutree2-c1-Z1_overDispersed_TPMmatrix.txt===== * This cluster is likely to be excitatory neurons (GAD1-/GAD2-/SLC6A1-) * 1st layer cut by cutree2, 2nd layer cut by either cutree or kmeans [[Image:20chip.TPMmatrix_0.2-0.7-Z1_cutree2-c1-Z1_overDispersed_TPMmatrix_cutree.jpg|600px]] [[Image:20chip.TPMmatrix_0.2-0.7-Z1_cutree2-c1-Z1_overDispersed_TPMmatrix_kmean2-5.jpg|600px]] =====20chip.TPMmatrix_0.2-0.7-Z1_kmean2-c1-Z1_overDispersed_TPMmatrix.txt===== * This cluster is likely to be excitatory neurons (GAD1-/GAD2-/SLC6A1-) * 1st layer cut by kmean2, 2nd layer cut by either cutree or kmeans [[Image:20chip.TPMmatrix_0.2-0.7-Z1_kmean2-c1-Z1_overDispersed_TPMmatrix_cutree.jpg|600px]] [[Image:20chip.TPMmatrix_0.2-0.7-Z1_kmean2-c1-Z1_overDispersed_TPMmatrix_kmean2-5.jpg|600px]] ==911 singles (13333 genes) passing the criteria of 40% relative hg19 percentage and 50% gene allowed missing== ===Overdispersed genes=== * 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=== * hierarchical clustering based on 20chip.TPMmatrix.newID_0.4filtered_0.5allowed_Z1_overDispersed_TPMmatrix.txt [[Image:HC_496.jpg|800px]] * CD figure [[Image:CD_496.jpg|800px]] * refined heatmap2 [[Image:HC-CD_496.jpg|800px]] ===K means clustering=== * K-mean clustering based on 20chip.TPMmatrix.newID_0.4filtered_0.5allowed_Z1_overDispersed_TPMmatrix.txt [[File:496_2c.jpg|400px]] [[File:496_3c.jpg|400px]] [[File:496_4c.jpg|400px]] [[File:496_5c.jpg|400px]]
Summary:
Please note that all contributions to ZhangLabWiki may be edited, altered, or removed by other contributors. If you do not want your writing to be edited mercilessly, then do not submit it here.
You are also promising us that you wrote this yourself, or copied it from a public domain or similar free resource (see
ZhangLabWiki:Copyrights
for details).
Do not submit copyrighted work without permission!
Cancel
Editing help
(opens in new window)
Template used on this page:
Template:Table
(
edit
)
Navigation menu
Personal tools
Not logged in
Talk
Contributions
Create account
Log in
Namespaces
Page
Discussion
English
Views
Read
Edit
View history
More
Search
Navigation
Main Page
Current events
Recent changes
Random page
Investigators
Matt Cai
Song Chen
Eric Chu
Dinh Diep
Elizabeth Duong
Shicheng Guo
Alan Fung
Daniel Jacobsen
Blue Lake
Huy Lam
Alice Li
Andrew Richards
Brandon Sos
Chris Wei
Yan Wu
Kun Zhang
Tools
What links here
Related changes
Special pages
Page information