讓你的單細胞數據動起來！iCellR二前情回顧：讓你的單細胞數據動起來！|iCellR(一)在上文我們已經基本完成了聚類分析和細胞在不同條件下的比例分析，下面繼續開始啦，讓你看看iCellR的強大！1. 每個cluster的平均表達量my.obj – clust.avg.exp(my.obj)head(my.obj@clust.avg)#gene cluster_1 cluster_2cluster_3 cluster_4 cluster_5 cluster_6#1 A1BG 0.072214723 0.092648973 0.08258609 0 0.027183115 0.072291636#2 A1BG.AS1 0.014380756 0.003280237 0.01817982 0 0.000000000 0.011545546#3 A1CF 0.000000000 0.000000000 0.00000000 0 0.000000000 0.000000000#4A2M 0.000000000 0.000000000 0.00000000 0 0.007004131 0.004672857#5A2M.AS1 0.003520828 0.039985296 0.00876364 0 0.056596203 0.018445562#6A2ML1 0.000000000 0.000000000 0.00000000 0 0.000000000 0.000000000# cluster_7cluster_8 cluster_9#1 0.09058946 0.04466827 0.027927923#2 0.00000000 0.01534541 0.005930566#3 0.00000000 0.00000000 0.000000000#4 0.00000000 0.00000000 0.003411938#5 0.00000000 0.00000000 0.000000000#6 0.00000000 0.00000000 0.000000000保存objectsave(my.obj錛 file = "my.obj.Robj")2. 尋找marker基因marker.genes – findMarkers(my.obj錛fold.change = 2錛padjval = 0.1) # 篩選條件dim(marker.genes)# [1] 1070 17head(marker.genes)#row baseMeanbaseSD AvExpInCluster AvExpInOtherClusters#LRRN3 LRRN3 0.01428477 0.1282046 0.055372430.003437002#LINC00176 LINC00176 0.06757573 0.2949763 0.214041510.028906516#FHIT FHIT 0.10195359 0.3885343 0.314049360.045957058#TSHZ2 TSHZ2 0.04831334 0.2628778 0.143009980.023311970#CCR7 CCR7 0.28132627 0.6847417 0.813864440.140728033#SCGB3A1 SCGB3A1 0.06319598 0.3554273 0.181305570.032013232#foldChange log2FoldChange pval padj clusters#LRRN316.110677 4.009945 1.707232e-06 2.847662e-031#LINC00176 7.404611 2.888424 4.189197e-16 7.117446e-131#FHIT6.833539 2.772633 1.576339e-19 2.681353e-161#TSHZ2 6.134616 2.616973 8.613622e-10 1.455702e-061#CCR75.783243 2.531879 1.994533e-42 3.400679e-391#SCGB3A1 5.663457 2.501683 2.578484e-07 4.313805e-041# genecluster_1 cluster_2 cluster_3 cluster_4 cluster_5#LRRN3 LRRN3 0.05537243 0.004102916 0.002190847 0 0.010902326#LINC00176 LINC00176 0.21404151 0.016772401 0.005203161 0 0.009293024#FHIT FHIT 0.31404936 0.008713243 0.022934924 0 0.035701186#TSHZ2 TSHZ2 0.14300998 0.008996236 0.009444180 0 0.000000000#CCR7 CCR7 0.81386444 0.075719109 0.034017494 0 0.021492756#SCGB3A1 SCGB3A1 0.18130557 0.039644151 0.001183264 0 0.000000000#cluster_6cluster_7 cluster_8 cluster_9#LRRN3 0.002087831 0.00000000 0.000000000 0.012113258#LINC00176 0.086762509 0.01198777 0.003501552 0.003560614#FHIT0.104189143 0.04144293 0.041064681 0.007218861#TSHZ2 0.065509372 0.01690584 0.002352707 0.015350123#CCR70.272580821 0.06523324 0.257130255 0.031304151#SCGB3A1 0.078878071 0.01198777 0.000000000 0.043410608# baseMean: average expression in all the cells 所有細胞的平均表達值# baseSD: Standard Deviation 標準偏差# AvExpInCluster: average expression in cluster number (see clusters) 該cluster中的平均表達值# AvExpInOtherClusters: average expression in all the other clusters在所有其他cluster的平均表達值# foldChange: AvExpInCluster/AvExpInOtherClusters表達差異倍數# log2FoldChange: log2(AvExpInCluster/AvExpInOtherClusters) 取表達差異倍數的log2的值# pval: P value# padj: Adjusted P value# clusters: marker for cluster number# gene: marker gene for the cluster cluster的marker gene# the rest are the average expression for each cluster單細胞分群後，怎麼找到Marker基因定義每一類群？3. 繪制基因表達情況基因MS4A1：tSNE 2D錛PCA 2D展示該基因在兩種降維方式下，每個組的平均表達分布圖。Box Plot和Bar plot展示了該基因在每個簇中的表達箱線圖和條形圖# tSNE 2DA – gene.plot(my.obj錛 gene = "MS4A1"錛plot.type = "scatterplot"錛interactive = F錛out.name = "scatter_plot")# PCA 2DB – gene.plot(my.obj錛 gene = "MS4A1"錛plot.type = "scatterplot"錛interactive = F錛out.name = "scatter_plot"錛plot.data.type = "pca")# Box PlotC – gene.plot(my.obj錛 gene = "MS4A1"錛box.to.test = 0錛box.pval = "sig.signs"錛col.by = "clusters"錛plot.type = "boxplot"錛interactive = F錛out.name = "box_plot")# Bar plot (to visualize fold changes)D – gene.plot(my.obj錛 gene = "MS4A1"錛col.by = "clusters"錛plot.type = "barplot"錛interactive = F錛out.name = "bar_plot")library(gridExtra)grid.arrange(A錛B錛C錛D) # 4個圖合並展示R語言 – 箱線圖（小提琴圖、抖動圖、區域散點圖）R語言 – 柱狀圖展示多個基因的plots:genelist = c("PPBP"錛"LYZ"錛"MS4A1"錛"GNLY"錛"LTB"錛"NKG7"錛"IFITM2"錛"CD14"錛"S100A9")###for(i in genelist){MyPlot – gene.plot(my.obj錛 gene = i錛interactive = F錛plot.data.type = "umap"錛cell.transparency = 1)NameCol=paste("PL"錛i錛sep="_")eval(call("-"錛 as.name(NameCol)錛 MyPlot))}###library(cowplot)filenames – ls(pattern="PL_")plot_grid(plotlist=mget(filenames[1:9]))熱圖分析:# 獲取top10高變基因MyGenes – top.markers(marker.genes錛 topde = 10錛 min.base.mean = 0.2錛filt.ambig = F)MyGenes – unique(MyGenes)# plotheatmap.gg.plot(my.obj錛 gene = MyGenes錛 interactive = T錛 out.name = "plot"錛 cluster.by = "clusters")# orheatmap.gg.plot(my.obj錛 gene = MyGenes錛 interactive = F錛 cluster.by = "clusters")R語言 – 熱圖美化4. 使用ImmGen進行細胞類型預測注意：ImmGen是小鼠基因組數據，此處的樣本數據是人的。對于157個ULI-RNA-Seq(Ultra-low-input RNA-seq)樣本錛使用的metadata是：https://github.com/rezakj/scSeqR/blob/dev/doc/uli_RNA_metadat.txt# 繪制cluster8中top40基因(平均表達值最小為0.2)在不同細胞的表達分布圖Cluster = 8MyGenes – top.markers(marker.genes錛 topde = 40錛 min.base.mean = 0.2錛 cluster = Cluster)# 畫圖cell.type.pred(immgen.data = "rna"錛 gene = MyGenes錛 plot.type = "point.plot")# andcell.type.pred(immgen.data = "uli.rna"錛 gene = MyGenes錛 plot.type = "point.plot"錛 top.cell.types = 50)# orcell.type.pred(immgen.data = "rna"錛 gene = MyGenes錛 plot.type = "heatmap")# andcell.type.pred(immgen.data = "uli.rna"錛 gene = MyGenes錛 plot.type = "heatmap")# And finally check the genes in the cells and find the common ones to predict# heatmap.gg.plot(my.obj錛 gene = MyGenes錛 interactive = F錛 cluster.by = "clusters")# 可以看出來cluster 8更像B細胞！# for tissue type prediction use this:#cell.type.pred(immgen.data = "mca"錛 gene = MyGenes錛 plot.type = "point.plot")cellassign：用于腫瘤微環境分析的單細胞注釋工具（9月Nature）5. Pathway analysis# cluster7的KEGG Pathway# pathways.kegg(my.obj錛 clust.num = 7)# this function is being improved and soon will be available 這個函數還要改進，

線上百家樂輸錢

很快就能用了Pathview包：整合表達譜數據可視化KEGG通路KEGG在線數據庫使用攻略對cluster進行QC：clust.stats.plot(my.obj錛 plot.type = "box.mito"錛 interactive = F)# 每個cluster中細胞的線粒體基因分布圖clust.stats.plot(my.obj錛 plot.type = "box.gene"錛 interactive = F)# 每個cluster中細胞的基因分布圖6. 差異表達分析diff.res – run.diff.exp(my.obj錛 de.by = "clusters"錛 cond.1 = c(1錛4)錛 cond.2 = c(2))# cluster之間的比較得到差異表達基因diff.res1 – as.data.frame(diff.res)diff.res1 – subset(diff.res1錛 padj 0.05) # 篩選padj 0.05的基因head(diff.res1)# baseMean1_4 2 foldChange log2FoldChange pval#AAK1 0.19554589 0.26338228 0.041792762 0.15867719-2.655833 8.497012e-33#ABHD14A0.09645732 0.12708519 0.027038379 0.21275791-2.232715 1.151865e-11#ABHD14B0.19132829 0.23177944 0.099644572 0.42991118-1.217889 3.163623e-09#ABLIM1 0.06901900 0.08749258 0.027148089 0.31029018-1.688310 1.076382e-06#AC013264.2 0.07383608 0.10584821 0.001279649 0.01208947-6.370105 1.291674e-19#AC092580.4 0.03730859 0.05112053 0.006003441 0.11743700-3.090041 5.048838e-07 padj#AAK1 1.294690e-28#ABHD14A1.708446e-07#ABHD14B4.636290e-05#ABLIM1 1.540087e-02#AC013264.2 1.950557e-15#AC092580.4 7.254675e-03# more examples 注意參數“de.by”設置的是不同差異比較方式diff.res – run.diff.exp(my.obj錛 de.by = "conditions"錛 cond.1 = c("WT")錛 cond.2 = c("KO"))#WT vs KOdiff.res – run.diff.exp(my.obj錛 de.by = "clusters"錛 cond.1 = c(1錛4)錛 cond.2 = c(2))# cluster 1 and 2 vs. 4diff.res – run.diff.exp(my.obj錛 de.by = "clustBase.condComp"錛 cond.1 = c("WT")錛 cond.2 = c("KO")錛 base.cond = 1)#cluster 1 WT vs cluster 1 KOdiff.res – run.diff.exp(my.obj錛 de.by = "condBase.clustComp"錛 cond.1 = c(1)錛 cond.2 = c(2)錛 base.cond = "WT") # cluster 1 vs cluster 2 but only the WT sample畫差異表達基因的火山圖和MA圖：# Volcano Plotvolcano.ma.plot(diff.res錛sig.value = "pval"錛sig.line = 0.05錛plot.type = "volcano"錛interactive = F)# MA Plotvolcano.ma.plot(diff.res錛sig.value = "pval"錛sig.line = 0.05錛plot.type = "ma"錛interactive = F)什麼？你做的差異基因方法不合適？7. 合並，重置，重命名和刪除cluster# 如果你想要合並cluster 2和cluster 3my.obj – change.clust(my.obj錛 change.clust = 3錛 to.clust = 2)# 重置為原來的分類my.obj – change.clust(my.obj錛 clust.reset = T)# 可以將cluster7編號重命名為細胞類型-"B Cell".my.obj – change.clust(my.obj錛 change.clust = 7錛 to.clust = "B Cell")# 刪除cluster1my.obj – clust.rm(my.obj錛 clust.to.rm = 1)# 進行tSNE對細胞重新定位my.obj – run.tsne(my.obj錛 clust.method = "gene.model"錛 gene.list = "my_model_genes.txt")# Use this for plotting as you make the changescluster.plot(my.obj錛 cell.size = 1錛 plot.type = "tsne"錛 cell.color = "black"錛 back.col = "white"錛 col.by = "clusters"錛 cell.transparency = 0.5錛 clust.dim = 2錛 interactive = F)Cell gating：可以框出指定的信息my.plot – gene.plot(my.obj錛 gene = "GNLY"錛plot.type = "scatterplot"錛clust.dim = 2錛interactive = F)cell.gating(my.obj錛 my.plot = my.plot)# or#my.plot – cluster.plot(my.obj錛#cell.size = 1錛#cell.transparency = 0.5錛#clust.dim = 2錛#interactive = F)下載細胞ID之後，用下面的命令重命名這些clustermy.obj – gate.to.clust(my.obj錛 my.gate = "cellGating.txt"錛 to.clust = 10)8. 偽時序分析注意參數“type”MyGenes – top.markers(marker.genes錛 topde = 50錛 min.base.mean = 0.2)MyGenes – unique(MyGenes)pseudotime.tree(my.obj錛 marker.genes = MyGenes錛 type = "unrooted"錛 clust.method = "complete")# orpseudotime.tree(my.obj錛 marker.genes = MyGenes錛 type = "classic"錛 clust.method = "complete")9. 應用monocle進行偽時序分析library(monocle)MyMTX – my.obj@main.dataGeneAnno – as.data.frame(row.names(MyMTX))colnames(GeneAnno) – "gene_short_name"row.names(GeneAnno) – GeneAnno$gene_short_namecell.cluster – (my.obj@best.clust)Ha – data.frame(do.call('rbind'錛 strsplit(as.character(row.names(cell.cluster))錛'_'錛fixed=TRUE)))[1]clusts – paste("cl."錛as.character(cell.cluster$clusters)錛sep="")cell.cluster – cbind(cell.cluster錛Ha錛clusts)colnames(cell.cluster) – c("Clusts"錛"iCellR.Conds"錛"iCellR.Clusts")Samp – new("AnnotatedDataFrame"錛 data = cell.cluster)Anno – new("AnnotatedDataFrame"錛 data = GeneAnno)my.monoc.obj – newCellDataSet(as.matrix(MyMTX)錛phenoData = Samp錛 featureData = Anno)## find disperesedgenesmy.monoc.obj – estimateSizeFactors(my.monoc.obj)my.monoc.obj – estimateDispersions(my.monoc.obj)disp_table – dispersionTable(my.monoc.obj)unsup_clustering_genes – subset(disp_table錛 mean_expression = 0.1)my.monoc.obj – setOrderingFilter(my.monoc.obj錛 unsup_clustering_genes$gene_id)# tSNEmy.monoc.obj – reduceDimension(my.monoc.obj錛 max_components = 2錛 num_dim = 10錛reduction_method = 'tSNE'錛 verbose = T)# clustermy.monoc.obj – clusterCells(my.monoc.obj錛 num_clusters = 10)## plot conditions and clusters based on iCellR analysisA – plot_cell_clusters(my.monoc.obj錛 1錛 2錛 color = "iCellR.Conds")B – plot_cell_clusters(my.monoc.obj錛 1錛 2錛 color = "iCellR.Clusts")## plot clusters based monocle analysisC – plot_cell_clusters(my.monoc.obj錛 1錛 2錛 color = "Cluster")# get marker genes from iCellR analysisMyGenes – top.markers(marker.genes錛 topde = 30錛 min.base.mean = 0.2)my.monoc.obj – setOrderingFilter(my.monoc.obj錛 MyGenes)my.monoc.obj – reduceDimension(my.monoc.obj錛 max_components = 2錛method = 'DDRTree')# order cellsmy.monoc.obj – orderCells(my.monoc.obj)# plot based on iCellR analysis and marker genes from iCellRD – plot_cell_trajectory(my.monoc.obj錛 color_by = "iCellR.Clusts")## heatmap genes from iCellRplot_pseudotime_heatmap(my.monoc.obj[MyGenes錛]錛cores = 1錛cluster_rows = F錛use_gene_short_name = T錛show_rownames = T)NBT|45種單細胞軌蹟推斷方法比較，110個實際數據集和229個合成數據集10. iCellR應用于scVDJ-seq數據# first prepare the files.# this function would filter the files錛 calculate clonotype frequencies and proportions and add conditions to the cell ids.my.vdj.1 – prep.vdj(vdj.data = "all_contig_annotations.csv"錛 cond.name = "WT")my.vdj.2 – prep.vdj(vdj.data = "all_contig_annotations.csv"錛 cond.name = "KO")my.vdj.3 – prep.vdj(vdj.data = "all_contig_annotations.csv"錛 cond.name = "Ctrl")# concatenate all the conditionsmy.vdj.data – rbind(my.vdj.1錛 my.vdj.2錛 my.vdj.3)# see head of the filehead(my.vdj.data)#raw_clonotype_id barcode is_cell contig_id#1 clonotype1 WT_AAACCTGAGCTAACTC-1True AAACCTGAGCTAACTC-1_contig_1#2 clonotype1 WT_AAACCTGAGCTAACTC-1True AAACCTGAGCTAACTC-1_contig_2#3 clonotype1 WT_AGTTGGTTCTCGCATC-1True AGTTGGTTCTCGCATC-1_contig_3#4 clonotype1 WT_TGACAACCAACTGCTA-1True TGACAACCAACTGCTA-1_contig_1#5 clonotype1 WT_TGTCCCAGTCAAACTC-1True TGTCCCAGTCAAACTC-1_contig_1#6 clonotype1 WT_TGTCCCAGTCAAACTC-1True TGTCCCAGTCAAACTC-1_contig_2#high_confidence length chainv_gene d_genej_gene c_gene full_length#1True693 TRA TRAV8-1 NoneTRAJ21 TRACTrue#2True744 TRBTRBV28TRBD1 TRBJ2-1TRBC2True#3True647 TRA TRAV8-1 NoneTRAJ21 TRACTrue#4True508 TRBTRBV28TRBD1 TRBJ2-1TRBC2True#5True660 TRA TRAV8-1 NoneTRAJ21 TRACTrue#6True770 TRBTRBV28TRBD1 TRBJ2-1TRBC2True#productive cdr3cdr3_nt#1 TrueCAVKDFNKFYFTGTGCCGTGAAAGACTTCAACAAATTTTACTTT#2 True CASSLFSGTGTNEQFF TGTGCCAGCAGTTTATTTTCCGGGACAGGGACGAATGAGCAGTTCTTC#3 TrueCAVKDFNKFYFTGTGCCGTGAAAGACTTCAACAAATTTTACTTT#4 True CASSLFSGTGTNEQFF TGTGCCAGCAGTTTATTTTCCGGGACAGGGACGAATGAGCAGTTCTTC#5 TrueCAVKDFNKFYFTGTGCCGTGAAAGACTTCAACAAATTTTACTTT#6 True CASSLFSGTGTNEQFF TGTGCCAGCAGTTTATTTTCCGGGACAGGGACGAATGAGCAGTTCTTC#reads umis raw_consensus_id my.raw_clonotype_id clonotype.Freq#112412 clonotype1_consensus_1clonotype1120#224004 clonotype1_consensus_2clonotype1120#310902 clonotype1_consensus_1clonotype1120#424554 clonotype1_consensus_2clonotype1120#513462 clonotype1_consensus_1clonotype1120#630738 clonotype1_consensus_2clonotype1120#proportion total.colonotype#1 0.04098361 1292#2 0.04098361 1292#3 0.04098361 1292#4 0.04098361 1292#5 0.04098361 1292#6 0.04098361 1292# add it to iCellR objectadd.vdj(my.obj錛 vdj.data = my.vdj.data)reference如果想嘗試一下，這裡有傳送門哦！https://github.com/rezakj/iCellRhttps://www.nature.com/articles/s41597-019-0351-8NGS系列文章包括NGS基礎、轉錄組分析 （Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析 （ChIP-seq基本分析流程）、單細胞測序分析 (重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step) – Limma差異分析、火山圖、功能富集）等內容。背景介紹腎髒是具有許多不同功能的高度復雜的器官，由幾個功能和解剖上不連續的部分組成。腎小球和腎小管是腎單位的重要組成部分。足細胞與腎小球內皮細胞一起合成了腎小球基底膜，

央視中秋晚會陣容

它是最終的過濾屏障，防止蛋白質損失到尿液中。頂葉上皮細胞（Parietal epithelial cells，PECs）是另一種常見的腎小球細胞類型，可能導致腎小球硬化、新月和假新月形成。近端小管（proximal tubule，PT）通過控制Na+ – H+和HCO3-的轉運在調節全身酸堿平衡中起著重要作用，而遠曲小管則更多地參與電解質的轉運。在先前的研究中，研究人員對腎髒不同組成部分進行了bulk RNA測序（RNA-seq最強綜述名詞解釋&思維導圖|關于RNA-seq，你想知道的都在這（續）），

金多電子遊戲場

為理解不同片段的轉錄組提供參考。然而，bulk RNA測序不能反映單細胞水平的轉錄組，只能反映總體平均RNA表達（自從用了這個神器，

小山村考出12位博士28位碩士

大規模RNA-seq數據挖掘我也可以）。正常人腎髒的全面細胞解剖結構對于解決腎髒疾病和腎癌的細胞起源至關重要。一些腎髒疾病可能是細胞類型特異性的，尤其是腎小管細胞。為了研究人腎髒的分類和轉錄組信息，作者迅速獲得了腎髒的單細胞懸液並進行了單細胞RNA測序（scRNA-seq）（重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述））。作者介紹了來自三個人類供體腎髒的23錛366個高質量細胞的scRNA-seq數據，並將正常人腎細胞劃分為10個clusters。其中，近端腎小管（PT）細胞被分為三個亞型，而集合導管細胞被分為兩個亞型。總體而言，該數據為腎細胞生物學和腎髒疾病的研究提供了可靠的參考。下面我們按照作者的分析思路復現該文章的部分內容：首先，從GSE131685中下載數據：裡面的文件名要分別改為“barcodes.tsv”、“genes.tsv”和“matrix.mtx”，在Read10X（Hemberg-lab單細胞轉錄組數據分析（七）- 導入10X和SmartSeq2數據Tabula Muris）時才不會報錯。。。Kidney data loadinglibrary(devtools)install_github("immunogenomics/harmony")library(Seurat)library(magrittr)library(harmony)library(dplyr)#Kidney data loading 並構建seurat objectK1.data – Read10X(data.dir = "/Users/zhanghu1992/Documents/GSE131685_RAW/kidney1/")K1 – CreateSeuratObject(counts = K1.data錛 project = "kidney1"錛 min.cells = 8錛 min.features = 200)K2.data – Read10X(data.dir = "/Users/zhanghu1992/Documents/GSE131685_RAW/kidney2/")K2 – CreateSeuratObject(counts = K2.data錛 project = "kidney2"錛 min.cells = 6錛 min.features = 200)K3.data – Read10X(data.dir = "/Users/zhanghu1992/Documents/GSE131685_RAW/kidney3/")K3 – CreateSeuratObject(counts = K3.data錛 project = "kidney3"錛 min.cells = 10錛 min.features = 200)kid – merge(x = K1錛 y = list(K2錛 K3)) #讀取文件並用merge函數進行合並插一句嘴，我們來看一下華盛頓大學PhD jared.andrews對merge函數的解釋：注意老鐵說的“Seurat’s integration method is quite heavy handed in my experience錛so if you decide to go the integration route錛I’d recommend using the SeuratWrapper around the fastMNN”（單細胞分析Seurat使用相關的10個問題答疑精選！）QC# quality controlkid[["percent.mt"]] – PercentageFeatureSet(kid錛 pattern = "^MT-") #提取有關線粒體的基因VlnPlot(kid錛 features = c("nFeature_RNA"錛 "nCount_RNA"錛 "percent.mt")錛 ncol = 3) #由圖可以看出分布還可以plot1 – FeatureScatter(kid錛 feature1 = "nCount_RNA"錛 feature2 = "percent.mt")plot2 – FeatureScatter(kid錛 feature1 = "nCount_RNA"錛 feature2 = "nFeature_RNA")CombinePlots(plots = list(plot1錛 plot2))kid – subset(kid錛 subset = nFeature_RNA 200 & nFeature_RNA 2500 & percent.mt 30) #篩選條件kid – NormalizeData(kid錛 normalization.method = "LogNormalize"錛 scale.factor = 10000)kid – NormalizeData(kid) #標準化kid – FindVariableFeatures(kid錛 selection.method = "vst"錛 nfeatures = 2000) #查找高變基因top10 – head(VariableFeatures(kid)錛 10)plot1 – VariableFeaturePlot(kid)plot2 – LabelPoints(plot = plot1錛 points = top10錛 repel = TRUE)CombinePlots(plots = list(plot1錛 plot2))# 計算細胞週期s.genes -cc.genes$s.genesg2m.genes-cc.genes$g2m.geneskid – CellCycleScoring(kid錛 s.features = s.genes錛 g2m.features = g2m.genes錛 set.ident = TRUE)all.genes – rownames(kid)kid – ScaleData(kid錛 vars.to.regress = c("S.Score"錛 "G2M.Score")錛 features = all.genes)這裡我想叨叨幾句，據我看到的文獻，多數是在進行降維後將細胞週期方面對分群的影響作為一個單獨模塊去敘述，作者在先期不管細胞週期對聚類是否有影響的情況下就對細胞週期相關基因進行去除也是比較明智的，因為作者並不想讓該因素混雜其中影響分群（如何火眼金睛鑑定那些單細胞轉錄組中的混雜因素）。#當然我們還是要看是否細胞週期真的有影響，

第二波藝fun券申請

感興趣的小伙伴可以看一下，確實是有一定影響的！#kid – ScaleData(kid錛 features = rownames(kid))#kid- RunPCA(kid 錛 features = c(s.genes錛 g2m.genes))#DimPlot(kid)利用harmony算法去除批次效應並細胞分類#Eliminate batch effects with harmony and cell classificationkid – RunPCA(kid錛 pc.genes = kid@var.genes錛 npcs = 20錛 verbose = FALSE)options(repr.plot.height = 2.5錛 repr.plot.width = 6)kid – kid %%RunHarmony("orig.ident"錛 plot_convergence = TRUE)　#等候時間較長，請溜達溜達吧harmony_embeddings – Embeddings(kid錛 'harmony')harmony_embeddings[1:5錛 1:5]kid – kid %%RunUMAP(reduction = "harmony"錛 dims = 1:20) %%FindNeighbors(reduction = "harmony"錛 dims = 1:20) %%FindClusters(resolution = 0.25) %%identity()new.cluster.ids – c(0錛1錛 2錛 3錛 4錛 5錛 6錛 7錛8錛9錛10)names(new.cluster.ids) – levels(kid)kid – RenameIdents(kid錛 new.cluster.ids)“harmony”整合不同平台的單細胞數據之旅鑑定Marker基因#Calculating differentially expressed genes (DEGs) and Save rds filekid.markers – FindAllMarkers(kid錛 only.pos = TRUE錛 min.pct = 0.25錛 thresh.use = 0.25)#尋找高變基因write.table(kid.markers錛sep=" "錛file="/home/yuzhenyuan/Seurat/0.2_20.xls")saveRDS(kid錛file="/home/yuzhenyuan/kid/har/0.25_20.rds")單細胞分群後，怎麼找到Marker基因定義每一類群？可視化#Some visual figure generationDimPlot(kid錛 reduction = "umap"錛 group.by = "orig.ident"錛 pt.size = .1錛 split.by = 'orig.ident'DimPlot(kid錛 reduction = "umap"錛 group.by = "Phase"錛 pt.size = .1)　#按照細胞週期進行劃分DimPlot(kid錛 reduction = "umap"錛 label = TRUE錛 pt.size = .1)　#注意作者在用同樣參數設置後分為10個clusters，其實無關緊要，

麻將怎麼拿牌

都需要通過marker重新貼現。根據作者提供的marker對細胞亞群進行貼現，如下圖所示：其實部分marker並不是特異性marker，所以在進行區分的時候一定要好好甄別。與以下原文圖基本相同，個人感覺tSNE是不是也有什麼隨機種子的東東，

央視中秋晚會陣容

感覺總會略有不同：DoHeatmap(kid錛 features = c("SLC13A3"錛"SLC34A1"錛"GPX3"錛"DCXR"錛"SLC17A3"錛"SLC22A8"錛"SLC22A7"錛"GNLY"錛"NKG7"錛"CD3D"錛"CD3E"錛"LYZ"錛"CD14"錛"KRT8"錛"KRT18"錛"CD24"錛"VCAM1"錛"UMOD"錛"DEFB1"錛"CLDN8"錛"AQP2"錛"CD79A"錛"CD79B"錛"ATP6V1G3"錛"ATP6V0D2"錛"TMEM213"))　#　繪制部分基因熱圖R語言 – 熱圖美化VlnPlot(kid錛 pt.size =0錛 idents= c(1錛2錛3)錛 features = c("GPX3"錛 "DCXR"錛"SLC13A3"錛"SLC34A1"錛"SLC22A8"錛"SLC22A7"))VlnPlot(kid錛 idents= c(8錛10)錛 features = c("AQP2"錛 "ATP6V1B1"錛"ATP6V0D2"錛"ATP6V1G3"))R語言 – 箱線圖（小提琴圖、抖動圖、區域散點圖）可視化之為什麼要使用箱線圖？tSNE plot##tSNE Plotkid -RunTSNE(kid錛 reduction = "harmony"錛 dims = 1:20)TSNEPlot(kid錛 do.label = T錛 label = TRUE錛 do.return = T錛 pt.size = 1)TSNEPlot(kid錛 do.return = T錛 pt.size = 1錛 group.by = "orig.ident"錛 split.by = 'orig.ident')TSNEPlot(kid錛 do.return = T錛 pt.size = 1錛 group.by = "Phase")Hemberg-lab單細胞轉錄組數據分析（十二）- Scater單細胞表達譜tSNE可視化與前面的圖是相同的。提取PT cells進行後續分析#Select a subset of PT cells（近端小管）PT – SubsetData(kid錛 ident.use = c(0錛1錛2)錛 subset.raw = T)，