‘壹’ RNAseq教程(2.2)
1.Mole 1 - Introction to RNA sequencing
2.Mole 2 - RNA-seq Alignment and Visualization
3.Mole 3 - Expression and Differential Expression
4.Mole 4 - Isoform Discovery and Alternative Expression
5.Mole 5 - De novo transcript reconstruction
6.Mole 6 - Functional Annotation of Transcripts
用HISAT2比对基因组和转录组。
首先，为对齐结果创建适当的输出目录
HISAT2的输出是每个数据集的SAM/BAM文件。
参考HISAT2帮助手册获得更多说明：
HISAT2基本用法:
额外参数如下:
注意:在上面的对齐中，我们将每个库视为一个独立的数据集。如果你有一个库的多个数据，你可以在一个HISAT2命令中将它们对齐在一起。要组合多个reads文件，您需要为'-1'输入参数提供所有read1文件作为逗号分隔的列表，然后为'-2'输入参数提供所有read2文件作为逗号分隔的列表，(其中两个列表的顺序相同):还可以使用samtool
将HISAT2 sam文件转换为bam文件，并按对齐位置排序
将所有UHR数据和所有HBR数据合并成一个BAM文件。注意:这可以通过几种方式来完成，比如“samtools merge”，“bamtools merge”，或者使用picard-tools(见下文)。我们选择第三种方法是因为它在合并bam头信息方面做得最好。注意:sambamba也保留头部信息。
计算对齐(BAM)文件，确保所有文件都成功创建(总共应该有8个)
任务:对额外的数据集执行一些比对。用你在上面学到的技巧来对齐阅读。尝试使用HISAT2。还要练习将SAM转换为BAM文件，并合并BAM文件。
在练习3中创建的名为“practice”的单独工作目录中进行分析。
.sam和.bam文件之间有什么区别?
sam文件是一个纯文本序列比对映射文件。bam文件是相同信息的二进制压缩版本。
如果您像上面所做的那样对结果BAM文件进行排序，那么结果是否按read名称排序?还是position?
按照position排序
可以查看BAM文件的哪些列以确定排序的样式?
第一、第三和第四列包含reads名称、染色体和位置。
可以使用什么命令仅查看BAM头?
‘贰’ RNA-seq原理
测序技术发展：
1977Sanger测序--1996焦磷酸测序--2003cmPCR--2003ZMW---2012纳米孔测序
RNA-seq的一些技术限制，测序误差主要由生物学误差(生物学重复，比如取30只小鼠采样)和技术性误差(技术性重复，比如对1只小鼠采样3次)造成，如果想要得到的数据为无偏的，那么生物学重复最重要，因为生物个体代表着样本，而技术手段只会造成不可控干扰。总的来说，只做技术性重复的实验结果偏差最大，技术性重复+生物学重复的实验结果偏差也可能较大，除非生物学重复远大于技术性重复(因为当生物学重复次数不足时，技术性重复能扩大样本单一的影响)，无论如何，多做生物学重复，这有助于你的结论被其他人复现。
原理详解：
A 为了保证细胞在标记的过程中是单独分开的，10X开发了微流体设备(microfuidic device)进行预处理，设备有三个上样孔，分别加入你的1.样本细胞悬液(Sample) 2.凝胶小球(Beads) 3.分离液(Oil)，下图为具体设备的示意图。
当我们把样本细胞悬液加入设备时，每一个细胞会与凝胶小球单独结合，然后被分离液包裹，形成一个油包水的密闭小液滴(droplet)。进一步地，细胞和凝胶小球相遇不久后会裂解，释放出里面的各种物质，RNA(mRNA、tRNA、rRNA)，蛋白质，脂质，DNA等。实际上Beads上联接了不同的接头，其中有一个接头包含ploy(dT)序列，在细胞裂解后释放的核酸中，只有mRNA带有polyA tail，于是Beads的poly(dT)接头就可以从众多的裂解产物里捕获到mRNA(实际上drop-seq采用3'端测序，就是为了检测polyA tail)。
Master Mix中带有反转录试剂，当mRNA被捕获后，就可以从它的3‘端开始作为模板，进行反转录出cDNA的第一条链，这第一条链就沿着poly(dT)序列延申，长在了beads上，形成了图一7中的STAMPs，接着我们把反转录出来的cDNA序列洗脱，以cDNA的第一条链为模板，进行PCR，合成cDNA的第二条链，然后就是我们熟悉的cDNA扩增以及illumina测序。
如何确定测序序列来自哪个细胞？single cell的RNA-seq和bulk的RNA-seq的最大区别是什么？是barcode，或者说是cell barcode(实际上DNA自带barcode，cell barcode是人为控制的)。每一种single cell的beads上都有着相同的cell barcode(beads与beads间的cell barcode是不同的)，假设每个beads只捕获一个cell，那么则每个cell都被cell barcode 单独标记了。
如何保证每个beads只捕获一个cell？第一是控制cell和beads的流速，第二是beads的数目远远超过cell的数目，即绝大多数的beads都是空的，只有少数的才捕获到了cell。但是还是有个别的droplet里面会两个或者更多的细胞，这就需要进行质控(QualityControl)。
接下来可以参照10X Genomics的说明书详解single cell RNA-seq的barcode。
实际上beads上一开始只接了Read1、Barcode、Poly(dT)。
名词解释：
Poly(dT): 用来和mRNA的polyA结合，捕获mRNA
UMI: 用来标记不同的PCR产物(用于count计数)。为了减少由于复制引起的误差(重复抽样导致重复计数)，人们在一些单细胞测序的步骤中增加了UMI(unique molecular identifiers)，UMIs 是由 4-10 个随机核苷酸组成的序列，在 mRNA 反转录后，进入到文库中，每一个 mRNA，随机连上一个 UMI，因此可以计数不同的 UMI，最终计数 mRNA 的数量。
10X Barcode: 用来标记不同的single cell
Sample Index: 用来标记不同的sample
P5和P7: 用来进行illumina的桥式PCR测序
Truseq Read 1、2: 用来进行连接beads，cDNA的PCR扩增和加P7接头
在这些序列中，P5、P7、Truseq Read 1、2 的序列是已知的。
其他的序列是怎么一步一步添加上去的？
具体步骤：
利用Poly(dT)来捕获mRNA，在mRNA的5'端插入TSO(Template Switch Oligo模板切换低聚糖)引物，然后从mRNA的polyA开始反转录，直至mRNA的DNA序列被转录完成，然后在beads序列的3'端插入CCC，再对mRNA的TSO进行反转录，至此完成了cDNA的第一条链(序列顺序和mRNA逆序)。上述步骤很重要，因为中间cDNA的序列我们是不知道的(仪器测序长度有限)，如果不加上这个接头，就没有办法设计引物来合成cDNA的第二条链。
将mRNA溶解，对cDNA的第一条链加入UMI引物，以cDNA的第一条链为模板合成cDNA的第二条链。最后使用PCR(聚合酶链式反应)对cDNA(拷贝DNA)进行扩增(为了富集)。
PCR原理
因为II代测序(NGS)的illumina测序不能测很长的seq，约为200-700bp，所以不能测得mRNA全长，因此需要进一步把合成的cDNA利用酶打断到illumina能测的长度(长度有些随机，比如300bp的cDNA能通过头尾150bp完整测序，但700bp的cDNA只能通过头尾150bp测序+参考基因组推断出来)。然后在cDNA的3'端插入Truseq Read2引物(和Truseq Read1引物匹配为头尾，中间序列就是reads)、P5、P7。
最后的测序数据(reads)从Truseq Read1后的10X Barcode开始，一直到Truseq Read2为止。
PCR扩增是对cDNA单链进行复制，后面的桥式PCR是对完整的样本进行复制(增加数据深度)，总的来说各个cDNA呈均匀分布，然后进行抽样。
RNA-seq plications有PCR plication(最主要)、cluster plication、optical plication。
实际上仪器会对核苷酸进行染色，然后判断颜色确定ATCG碱基，因此有很多原因会导致机器误判，和后续QC有关。
1.某些核苷酸对颜色附着不明显
2.大片区域颜色相同(相同类型核苷酸)，而其中仅有几个颜色不同的点(不同类型的核苷酸)
‘叁’ 单细胞测序（sc-RNA seq）分析：Seurat4.0系列教程、高级分析、文章复现
时代的洪流奔涌而至，单细胞技术也从旧时王谢堂前燕，飞入寻常百姓家。雪崩的时候，没有一片雪花是无辜的，你我也从素不相识，到被一起卷入单细胞天地。R语言和Seurat已以势如破竹之势进入4.0时代，天问一号探测器已进入火星乌托邦平原了，你还不会单细胞吗？那么为了不被时代抛弃的太远，跟着我们一起通过学习seurat系列教程入门单细胞吧。
首先我们学习经典的标准流程，这个教程我当初入门时候曾经花1000购买过，也曾花6000购买过，不同的单细胞培训班，买的是一样的教程。现在免费送给你，别说话，开始学吧！
Seurat4.0系列教程1：标准流程 - (jianshu.com)
Seurat4.0系列教程2：多模式数据联合分析 - (jianshu.com)
Seurat4.0系列教程3：合并两个样品的10x数据集 - (jianshu.com)
Seurat4.0系列教程4：整合分析 - (jianshu.com)
Seurat4.0系列教程5：交互技巧 - (jianshu.com)
Seurat4.0系列教程6：常用命令 - (jianshu.com)
Seurat4.0系列教程7：数据可视化方法 - (jianshu.com)
Seurat4.0系列教程8：细胞周期评分和回归分析 - (jianshu.com)
Seurat4.0系列教程9：差异表达检测 - (jianshu.com)
Seurat4.0系列教程10：降维 - (jianshu.com)
Seurat4.0系列教程11：使用sctransform - (jianshu.com)
Seurat4.0系列教程12：大数据集整合的方法 - (jianshu.com)
Seurat4.0系列教程13：使用互惠 PCA （RPCA） 快速整合数据 - (jianshu.com)
Seurat4.0系列教程14：整合scRNA-seq and scATAC-seq数据 - (jianshu.com)
Seurat4.0系列教程15：映射和注释查询数据集 - (jianshu.com)
Seurat4.0系列教程16：多模式参考映射注释细胞 - (jianshu.com)
Seurat4.0系列教程17：Mixscape - (jianshu.com)
Seurat4.0系列教程18：Weighted Nearest Neighbor Analysis（WNN)） - (jianshu.com)
Seurat4.0系列教程19：多线程并行策略 - (jianshu.com)
Seurat4.0系列教程20：单细胞对象的格式转换 - (jianshu.com)
Seurat4.0系列教程21： 结合Cell Hashing分析双细胞 - (jianshu.com)
Seurat4.0系列教程22：空间转录组的分析、可视化与整合 - (jianshu.com)
Seurat4.0系列教程告一段落，但这决不是终点。这个系列教程只是给大家打开一扇窗，知道Seurat4.0有这些功能可用，少走弯路，起到一个抛转引玉的作用，后续还要自己深入研究。不要像我当初入门单细胞之时，为了找到整合方法耗费大量时间及不必要的金钱。君子生非异也，善假于物也。学，然后知不足。加油吧！少年！用好单细胞测序及分析这个技术，为人类健康研究做出自己的贡献！不足之处，恳请批评指正！------以此纪念粉丝突破1000。
CellChat三部曲1：使用CellChat对单个数据集进行细胞间通讯分析 - (jianshu.com)
CellChat 三部曲2：使用CellChat 对多个数据集细胞通讯进行比较分析 - (jianshu.com)
CellChat 三部曲3：具有不同细胞类型成分的多个数据集的细胞通讯比较分析 - (jianshu.com)
RNAvelocity系列教程1：RNA速率简介及scVelo安装 - (jianshu.com)
RNAvelocity系列教程2：使用 Seurat 和 scVelo 估计 RNA 速率 - (jianshu.com)
RNAvelocity系列教程3：使用Seurat和velocyto估算RNA速率 - (jianshu.com)
RNAvelocity系列教程4：velocyto.R的使用教程 - (jianshu.com)
RNAvelocity系列教程5：scVelo 的基本知识 - (jianshu.com)
RNAvelocity系列教程6：# scVelo用于RNA 速率基本流程 - (jianshu.com)
RNAvelocity系列教程7：# scVelo实战-胰腺数据集 - (jianshu.com)
RNAvelocity系列教程8：# scVelo实战-齿状回数据集 - (jianshu.com)
RNAvelocity系列教程9：# scVelo应用-动力学模型 - (jianshu.com)
RNAvelocity系列教程10：# scVelo应用-微分动力学 - (jianshu.com)
使用fgsea进行单细胞转录组GSEA富集分析 - (jianshu.com)
Seurat24节气之21大雪---vars.to.regress = "percent.mt"排除线粒体基因 - (jianshu.com)
Seurat24节气之2雨水---resolution参数 - (jianshu.com)
单细胞转录组基础分析五：细胞再聚类 - (jianshu.com)
单细胞测序文章图表复现01—读入文件的表达矩阵构建Seurat对象 - (jianshu.com)
单细胞测序文章图表复现02—Seurat标准流程之聚类分群 - (jianshu.com)
单细胞测序文章图表复现03—区分免疫细胞和肿瘤细胞 - (jianshu.com)
单细胞测序文章图表复现04—聚类分群的resolution参数问题 - (jianshu.com)
CNS图表复现11—RNA-seq数据可不只是表达量矩阵结果 - (jianshu.com)
‘肆’ 单细胞综述之整合分析
文章发表于nature review genetics： Integrative single- cell analysis ，作者是Tim Stuart与 Rahul Satija 。做过单细胞分析的对他们应该不陌生。
scRNA-seq技术的发展契合了研究个体细胞表观遗传、空间研究、蛋白质组与谱系信息的方法需要，这为研究多类型数据的综合方法提出了独特的机遇与挑战。综合分析可以发现细胞之间的模式关系，获取细胞的整体状态信息，产生涵盖不同样本与不同研究手段的数据集。该文重点讨论了单细胞基因表达数据与其他类型的单细胞分析方法的整合。
多模态（Multimodal）数据 ：多种类型数据的组合，如RNA与蛋白质数据组合，是一种多维度数据，类似多组学。
单模态 ：单个类型数据
Pseudotime ：拟时分析
联合聚类（Joint-clustering） ：通过联合不同类型数据对细胞进行分组。
典型相关分析（CCA） ： 利用综合变量对之间的相关关系来反映两组指标之间的整体相关性的多元统计分析方法。
动态时间规整（Dynamic time warping） ：一种局部拉伸或压缩两个一维矢量以校正一个矢量相对于另一个矢量的滞后的方法。
MNNs ：标准化基因表达空间中最临近的细胞。聚类用校正批次效应。
梯度推进（Gradient boosting） ：一种预测模型算法。
随着分子生物学、微流控与纳米技术的发展，催生了许多类型的单细胞测序技术。过去的方法集中在单模态测量上，如DNA序列、RNA表达量和 染色质可及性 上。虽然这些技术促进了我们对细胞多样性与发育景观的理解，但是它们并不能很好地解析单细胞内分子间互作关系。而这些互作关系是深入探索细胞状态的关键。此外，随着可用数据集规模的快速增长，迫切需要用于标准化与联合分析且考量到批次效应与个体差异的计算方法。
scRNA-seq是应用最为广泛的单细胞测序技术之一。而后出现了一系列互补技术如单细胞基因组、表观基因组和蛋白质组分析技术，涵盖了单细胞基因组测序（ Vitak, S. A. et al., 2017 ; Navin, N. et al., 2011 ）、染色质可及性（ Pott, S., 2017 ; Corces, M. R. et al., 2016 ; Buenrostro, J. D. et al., 2015 ; Cusanovich, D. A. et al., 2015 ; Lake, B. B. et al., 2018 ）、DNA甲基化（ Luo, C. et al., 2017 ; Smallwood, S. A. et al., 2014 ; Guo, H. et al., 2013 ; Mulqueen, R. M. et al., 2018 ）、膜蛋白（ Stoeckius, M. et al., 2017 ; Peterson, V. M. et al., 2017 ）、小RNA（ Faridani, O. R. et al., 2016 ）、组蛋白修饰（ Gomez, D. te al., 2013 ; Rotem, A. et al., 2015 ）和染色体构象（ Ramani, V. et al., 2017 ; Nagano, T. et al., 2013 ）等技术。目前已开发出研究单细胞空间结构和谱系信息的方法（ Frieda, K. L. et al., 2017 ; Shah, S. et al., 2016 ）。
单细胞多模态综合分析方法示意
单模态与多模态分析方法汇总
CEL-seq ：线性扩增测序法
CITE- seq ：膜蛋白丰度与基因表达水平测定
G＆T-seq ：基因组转录组测序
LINNAEUS ：谱系追踪
MARS-seq ：大规模平行单细胞RNA测序
MEMOIR ：谱系与空间结构测定
MERFISH ：主要是细胞间结构测定
osmFISH ：环状单分子荧光原位杂交，空间结构测定
REAP- seq ：膜蛋白丰度与基因表达水平测定
scATAC-seq ：单细胞空间结构测定
scBS-seq ：单细胞甲基化测序
scChIP-seq ：单细胞ChIP-seq
scGESTALT ：结合CRISPR-cas9的谱系追踪弄方法
scHi-C-seq ：测定染色体组装
sciATAC-seq ：结合index转座酶的scATAC-seq
sci-CAR ：利用index联合分析mRNA和染色质可及性谱
sci-MET ：利用index分析单细胞甲基化水平
sci-RNA-seq ：结合index的scRNA-seq
SCI-seq ：单细胞组合标记测序，检测CNV
scM＆T-seq ：单细胞甲基化组和转录组测序，可研究未知的DNA甲基化与基因表达之间的关系
scNOMe- seq ：核小体占位与甲基化组测序
scRRBS ：单细胞限制性代表区域甲基化测序
scTHS- seq ：单细胞转座体超敏性位点测序
seqFISH ：内含子序贯荧光原位杂交，扩展观测到基因数量
snmC-seq ：单核甲基胞嘧啶测序
SNS ：单核测序
SPLiT-seq ：丐版scRNA-seq
STARmap ：原位单细胞测序
理想的实验流程应当全面洞悉细胞的所有方面，包括分子状态、空间构象、胞外环境互作的全部过程。尽管当下技术手段无法做到，但多模态技术与综合计算方法可以是我们离该目标越来越近。文章希望提出整合单细胞转录组学、基因组学、表观组学与蛋白组学的数据统一分析方法，重点在结合其他数据类型分析scRNA-seq数据，尤其是整合来自于同一细胞的不同类型数据。
文章分为四大块，首先探讨了多模态单细胞分析方法，其次研究了不同实验不同数据整合分析，然后讨论了单细胞空间测序数据整合分析方法，最后给出了整合分析方法的前景与必要性。
最初的单细胞分析方法主要关注细胞某状态下的某类分子水平。而现在更引人瞩目的是同时分析单细胞内多种分子以建立更全面的单细胞分子视图。通常这些方法是将scRNA-seq数据与其它分析手段的结合，目前主要有四种策略从单细胞中得到多模态数据：
严格来说这种方法算单模态。
一些scRNA-seq workflow采用流式分选细胞，随后进行scRNA-seq（MARS-seq/Smart-seq/2），这样可以同时获得单细胞与对应的荧光信号，将荧光所表示的蛋白质水平与转录组在同一细胞中关联（ Ramsköld, D. et al., 2012 ; Jaitin, D. A. et al., 2014 ; Picelli, S. et al., 2013 ）。早期研究（ Hayashi, T. et al., 2010 ）利用FACS结合半定量RT-PCR（作者称之为FBSC‐PCR），结合scRNA-seq，明确了细胞表面marker可以区分细胞类型与状态( Wilson, N. K. et al., 2015 ;该文结合了Smart-seq2)，( Paul, F. et al., 2015 ;该文结合了MARS-seq)和鉴定稀有细胞的思路。 Paul, F. et al., 2015 与 Nestorowa, S. et al., 2016 利用该workflow研究发现了小鼠造血祖细胞由转录组定义不同细胞簇的免疫表型， Wilson, N. K. et al., 2015 则分离了小鼠HSCs，鉴定细胞维持干性相关的表面marker。但是囿于荧光光谱的重叠现象，利用该法测到的每个细胞的参数范围有限。
针对荧光无法分选的部分，FACS显然是不合适的，尤其是需要同时测得单细胞基因组与胞内蛋白的scRNA-seq实验。此时需要物理分离或通过不同tag筛选出不同组分。
G&T-seq通过加入oligo(dT)特异性分离mRNA同时保留基因组DNA从而实现了基因组转录组平行测序（ Macaulay, I. C. et al., 2015 ）DR-seq通过则通过加入barcode特异扩增cDNA序列实现基因组转录组平行测序（ Dey, S. S. et al., 2015 ）。这使得单细胞基因表达水平与其对应基因型联系起来，深度揭示单细胞间DNA拷贝数变异与染色体重排对下游mRNA丰度的具体关联。这些方法适用于研究体细胞基因高度变异的肿瘤组织。
DNA甲基化与转录组水平结合研究是基于 Macaulay, I. C. et al., 2015 的G&T-seq和 Smallwood, S. A. et al., 2014 的scBS- seq技术发展的，同普通BSP一样，用亚硫酸氢钠处理DNA片段随后进行扩增，结合G&T-seq，可以分析同一细胞内的DNA甲基化模式和基因表达数据（ Angermueller, C. et al., 2016 ）。由于DNA甲基化存在不稳定性和异质性，因此若要研究DNA甲基化与基因表达间的关系，则必须将表观基因组变异与细胞间的异质性区别开来。
通过DNA甲基化与转录组关联分析，为启动子甲基化与基因表达间的负相关性提供深层次的证据。此外，利用barcode系统选择性标记基因组DNA与cDNA，结合index系统，可以对数千个单细胞进行染色质可及性与基因表达水平间的关联分析，同时鉴定出影响基因表达的顺式调控元件（ Cao, J. et al., 2018 ）。
关于胞内蛋白与mRNA关联研究，有两种思路可供借鉴。其一（ Darmanis, S. et al., 2016 ）是将FACS sort到的细胞裂解后分离裂解液，分别进行蛋白质与RNA定量。作者采用 PEA (邻近探针延伸分析) 检测蛋白并用RT-qPCR定量，采用qRT-PCR定量mRNA。该法可以同时检测82个mRNA/75个蛋白；其二（ Genshaft, A. S. et al. ）是将FACS sort到的细胞在微流控芯片中同时进行逆转录和PEA而不分离裂解液。该法可以同时检测96个mRNA/38个蛋白。这两种方法检测的蛋白与mRNA数量与质量均有限。
这些技术的出现表明若将可以细胞信息转化为有序的barcode，我们就可以在分析单细胞转录组时将这些信息同时获取。这种策略不仅适用于分析细胞的自然状态，也适用于大规模基因扰动研究。目前有Perturb-Seq（ Dixit, A. et al., 2016 ）和CRISPR-Seq（ Adamson, B. et al., 2016 ; Datlinger, P. et al., 2017 ; Jaitin, D. A. et al., 2016 ），他们将scRNA-seq与CRISPR-cas9结合进行遗传筛选，使得研究正向遗传学的大规模基因扰动试验成为可能。具体原理是给单个基因扰动和受到影响的细胞添加barcode，通过scRNA-seq能够鉴定出这两者，从而推断CRISPR靶向基因和由此产生的单个细胞的转录谱间的关系。目前应用在基因调控网络（ Dixit, A. et al., 2016 ）、未折叠蛋白反应（ Adamson, B. et al., 2016 ）、免疫细胞分化发育（ Datlinger, P. et al., 2017 ）和T细胞受体激活（ Jaitin, D. A. et al., 2016 ），非编码区调控元件（ Klann, T. S. et al., 2017 ）。此外，还可以结合CRISPR-dcas9系统，扩展到转录调控、表观遗传调控领域中（ Thakore, P. I. et al., 2016 ; Liu, X. S. et al., 2016 ; Hilton, I. B. et al., 2015 ; Konermann, S. et al., 2015 ; Gilbert, L. A. et al., 2017 ）,18年发展了同时靶向和敲除基因的技术（ Boettcher, M. et al., 2018 ）。
另一个应用是结合CRISPR-cas9的谱系追踪技术。单细胞谱系追踪是去年的大热方向之一，此处提到三种mRNA+lineage方法： scGESTALT 、 ScarTrace 、 LINNAEUS 。这三种方法各有不同，但大体是利用CRISPR-cas9连续切割结合到基因组上的barcode，细胞会用NHEJ来应对这种损伤。但NHEJ容易出错，从而在DNA序列中产生随机突变，这些突变通过细胞分裂进行遗传，结合scRNAseq利用这些突变作为复合barcode来构建组织或器官发育谱系。
另一种略有不同的方法是 MEMOIR ，它结合smFISH与CRISPR-cas9系统，可以同时检测细胞谱系与空间位置。
普通的scRNA-seq流程除了可以做转录本丰度外，还可以进行诸如体细胞突变、遗传变异、RNA isoform等分析。
关于体细胞突变目前已有研究（ Lodato, M. A. et al., 2015 ）,该文通过对人大脑的少量单细胞全基因组测序，分析了发生的细胞突变，构建了人大脑神经细胞谱系。作者发现突变大多发生在高转录活性相关位置，这表明可能可以通过scRNA-seq数据来分析神经细胞突变情况，根据转录状态重构神经细胞谱系。此外，分析scRNA-seq数据中的拷贝数变异，可以研究癌症非整倍体与异质性等情况（ Tirosh, I. et al., 2016 ; Fan, J. et al., 2018 ）。
单细胞分析也为理解DNA自然变异如何影响基因表达与细胞状态提供了新思路。有研究结合GWAS+scRNAseq，鉴定出了不同个体之间的eQTL（ Kang, H. M. et al., 2018 ）。
多模态测序策略正在催生与之相匹配的数据分析方法。多模数据集可以检测到细胞间的细微差异，而单模数据很可能无法做到这一点。由于scRNAseq数据存在dropout，故而它更容易忽略细胞间的细微差别；但与来自同一细胞的其他数据互补分析可以改善这一问题。例如，很难通过scRNA-seq数据区分不同的T细胞亚群，但联合膜蛋白分析则可以显着提高亚群分辨率（ Stoeckius, M. et al., 2017 ），同样，RNA+chromatin、RNA+methylation联合可能揭示单个细胞间的调控异质性，不再赘述。
单细胞多模态分析思路很可能受到bulk-seq多组学联合分析的启发（ Meng, C. et al., 2016 ）， Argelaguet 开发了一种名为MOFA( multi- omics factor analysis)的方法，该方法在多组学bulk-seq数据中效果良好，同时测试了单细胞DNA甲基化数据与RNA数据联合处理情况，效果也可以。这暗示适用于bulk-seq的多组学数据处理方式可能也适用于单细胞多模态数据。鉴于单细胞数据规模远超bulk-seq，多视图机器学习不失为一种重要的补充手段（ Colomé- Tatché, M. & Theis, F. J., 2018 ）。
单细胞多模态研究策略为解析细胞内不同组分间的关系提供了新方法。如CITE-seq和REAP-seq可以轻易鉴别出相关度较低的RNA-protein模块，表明此处存在活跃的转录后调节。还有一个很有意思的是通过测量剪接过的成熟RNA与未剪接RNA的相对丰度，可以建立RNA与蛋白的关联动态模型（ La Manno, G. et al., 2018 ）。
此外，还可以在不同类型数据间建立统计模型。前面提到的sci-CAR文章建立了染色质可及性与基因表达水平间的统计模型，通过染色质可及性数据估计细胞内基因表达水平（ Cao, J. et al., 2018 ），另一组研究人员建立了gRNA与基因表达水平间的线性回归模型，用以识别细胞应答的前后关系，重构转录网络（Perturb-Seq（ Dixit, A. et al., 2016 ））。通过这种手段可以研究目标物种复杂的调控网络。
前面主要讲了在同一测序实验同一批细胞进行的多模态数据整合，而不同测序实验数据整合分析才是亟需解决的关键问题。同bulk seq 数据一样，处理批次效应是综合分析不同实验室、不同workflow产出数据的首要问题（SVA包（ Leek, J. T. 2014 ））。然而目前bulk seq水平的处理方法无法处理单细胞数据（（ Haghverdi, L, et al., 2018 ，作者用MNN处理数据，该法在 mnnpy 中得到改进）; Butler, A, et al,. 2018 ）。目前最新方法利用 CCA / MNN 可以识别出两个数据集间共有的部分，判定细胞间共有的生物学状态,然后以这些相同状态的细胞为基准消除批次效应。
此处作者介绍了他自己在Seurat V2中开发的方法（ Satija, R, et al., 2015 ;），该法用 CCA 鉴别出不同数据集间相同的细胞类型且可以避免出现由批次效应或常规PCA造成的假阳性细胞类型；接下来采用动态时间规整算法校正数据集间细胞密度差异。这两步骤可以将细胞投影到一个低维空间，具有相同生物学状态的细胞相互接近且消除了不同数据集带来的影响。
另一种方法即mnnCorrect，最早用于计算机领域图形识别。该法寻找不同数据集间最接近的细胞，将之判定为潜在的状态相同细胞，随后利用成对MNNs距离计算一个批次参数（batch vector），用以校正原始表达矩阵（ Haghverdi, L., 2018 ）。
CCA/mnnCorrect在整合处理不同来源的scRNA-seq数据时表现良好。这将极大提升发现稀有细胞、微弱转录差异细胞及与之对应maker的能力（ Haghverdi, L, et al,.2018 ； Butler, A,et al,. 2018 ） 。这为建立一个统一的单细胞参考数据集提供了依据。在此基础上，scRNA-seq数据整合分析得到了快速发展（ Hie, B. L, et al., 2018 ; Barkas, N. et al., 2018 ; Park, J.-E., 2018 ; Korsunsky, I. et al., 2018 ; Stuart, T. et al., 2018 ; Welch, J. et al., 2018 ）。这种多数据集整合分析的应用远不止用于校正批次效应这么单一。它可以在单细胞尺度上深入比较细胞间的状态，发现细胞对环境及基因扰动的特异性响应，对不同疾病及不同治疗下的患者的测序数据进行标准化。
scRNA-seq数据整合分析还可以扩展至跨物种分析。 Karaiskos,N 比较了两种果蝇早期胚胎的空间基因表达模式，通过构建空间基因表达图谱，该研究系统比较了两个果蝇的同源基因表达谱，鉴定出了彼此间的进化波动。 Tosches 比较了爬行动物与哺乳动物脑细胞间的相关性。 Baron 分析了人与小鼠胰岛细胞scRNA-seq数据，鉴定出了二者间的保守亚群。 Alpert 开发出了cellAlign，在一维水平上比对了人与小鼠的拟时轨迹，发现人胚胎合子激活要比小鼠晚，小鼠中比人活跃的基因皆与蛋白合成相关。跨物种分析未来是光明的，但对于多物种整合分析而言，精确鉴定物种间同源基因是多物种整合分析至关重要的一步。
以细胞分类信息的形式串联不同的scRNA-seq数据集，或者借鉴到自己实验中，是优于合并数据集然后de novo聚类这种方法的。且随着 有参细胞图谱 的开发，这种方式将更加寻常。目前已开发对应方法： scmap- cell & scmap- cluster ，其中scmap-cell 用乘积量化（ proct quantization ）算法进行比对，而scmap-cluster则用于识别未知数据集中的cluster。
利用已有的注释数据集，目前开发的新方法采用 奇异值分解 、 线性判别分析 和 支持向量机 算法来对细胞进行分类。此外，随着引用数据集的大小、范围与深度越来越高，监督聚类在解析细胞类型方面要比无监督聚类强得多。通过以上这些方法，可以更精确地识别并解析细胞亚群。
satija已有相关文章研究： Comprehensive Integration of Single-Cell Data
这一部分讲的是将scRNA-seq数据与其它不同来源和类型数据诸如甲基化、染色质结构等整合分析的方法。
将scRNA-seq数据与其它类型、不同来源的单细胞数据整合分析是无法提取到数据间的共同特征的，因为它们不是一个类型的数据，需要不同的分析方法。这点在基于基因组的数据（如染色质可及性与甲基化数据）与基于基因的数据（如基因与蛋白表达数据）间整合分析尤为明显。但如果这些数据来自于同一类细胞群，由于存在着共同的生物学状态，此时可以联立分析以发现不同数据集类型间的对应关系。
MATCHER 是一种在一维水平上比较不同类型测序数据拟时轨迹的方法。简单来说就是比对不同类型测序数据的拟时轨迹，以确定这些数据集间的对应关系。这种方法可以识别不同数据集间的“等效细胞”而不需预先知道彼此间的对应关系。开发者用scM&T- seq（ Angermueller, C. et al., 2016 ）和scRNA-seq数据做了验证，准确预测了DNA甲基化与基因表达之间的关系。
其他sc-seq数据不同于scRNA-seq数据一样可以借助Marker解析细胞类型，因此可以利用scRNA-seq解析出的细胞信息为其他sc-seq数据分析做参考。有研究（ Lake, B. B. et al., 2018 ）对不同脑组织切片进行了单核RNAseq（snRNA-seq）与单细胞转座子超敏性位点测序（scTHS-seq），通过梯度推进算法利用单细胞基因表达谱指导了染色质可及性测序数据集的细胞分类：作者首先鉴别出snRNA-seq数据集与scTHS-seq数据集共有的细胞亚群，训练一个可以将基因表达与染色质可及性数据关联的模型；然后利用该模型去分类scTHS-seq中剩余未被分类的细胞。这种方法可以更细致地对大脑组织中的细胞进行分类。同样，可以整合scATAC-seq数据集来分析单细胞DNA甲基化或转座酶染色质可及性间的细胞分类。
目前正在开发的新方法有利用假定等价特征、或识别在所有类型数据中的假定相关共享特征来进行数据交叉模态分类。 Welch 开发了一种集成非负矩阵分解（iNMF）的方法，名为LIGER，可以跨模态整合数据。他们对同一类型 皮质细胞 分别进行了亚硫酸盐测序（snmC- seq）与scRNA-seq并对其进行了分类。他们假设基因体甲基化与其表达水平负相关从而整合了不同模态测序数据进行细胞分类。在seurat v3.0中，作者也引入了假定等价特征或关联特征进行多模态整合数据细胞分类的方法。这些方法优点如上所述，即可以利用scRNA-seq的细胞分类信息来指导scATAC-seq数据细胞分类，鉴别出染色质可及性与DNA甲基化的细胞特异模块。
组织中细胞的空间结构常反映出细胞间的功能差异与细胞命运和谱系的差异。不同基因表达引导细胞向不同方向分化，不同细胞精确排列形成不同组织。关键是单细胞实验通常在分析前细胞已被解离，组织原位信息无法保留，scRNA-seq得到的表达谱不能完全反应细胞空间信息。具有相似基因表达谱的细胞可能存在于不同的空间位置中，故而细胞分离过程中空间信息的缺失是很多单细胞实验的主要缺点。结合高分辨率基因表达谱与空间表达图谱 (spatial expression maps) 将细胞空间坐标与基因表达谱联系起来，可以解决这一问题。有两类方法：计算模型或者RNA原位定量，可以同时收集到细胞空间坐标与基因表达值。
‘伍’ RNA-seq 分析之我见（一）
先说下生物体内RNA的大致组成：
编码RNA:根据中心法则我们知道，DNA转录为mRNA，mRNA通过tRNA翻译为蛋白质，蛋白质行使生命功能,例如呼吸，运动，消化等等。人类只有2万左右个蛋白质编码基因，这些编码基因只占人类全基因组的2%左右。mRNA占细胞RNA总量的2%～5%， tRNA占细胞RNA总量的15%左右。
非编码RNA:有些DNA转录为RNA后，不继续编码蛋白质，这种RNA叫非编码RNA（ncRNA），包括microRNA,lncRNA,cirRNA，之前人们认为这些RNA是“垃圾”，但是近年来的研究证明，这些RNA对编码基因发挥着重要的调控作用，是当下研究的热点。
rRNA：核糖体RNA,占RNA总量的80%左右。
广义上说占总RNA95%左右的rRNA和tRNA也属于非编码RNA，但是一般研究中，使用的是它的狭义概念，即除去rRNA和tRNA之外的非编码RNA。
正常情况下，非编码RNA调控基因的转录翻译，这些都是有序进行的。
但是当处于异常条件下，或者由于自身衰老变异或者受到外部的刺激，比如细菌病毒的感染，射线照射等，这之后往往导致非编码RNA表达的变化，进而影响蛋白表达的变化，从而引起一系列的病理反应，最终导致疾病。
那么反过来，如果我们想了解某一疾病具体的发病机理，我们是不是可以提取某一疾病状态下组织或者细胞的总RNA，去分析它们和正常组表达的异同，我们有理由相信，这些差异表达的RNA分子，很可能与发病机制有关，研究这些差异分子，可以给我们对这一疾病的发病机制的研究提供重要线索，从而研发出更有效的诊断和治疗方法。
通过上面的分析，接下来面临的问题就是，我怎么分析某一疾病状态下组织或者细胞所有RNA的表达情况，一个一个分析，肯定不现实，而且可能还有很多未被发现但是很重要的分子。怎么办？只有一个办法，转录组测序，即RNA-Seq, 某一条件下所有转录出来的RNA碱基序列，我都给你测出来是什么。
那么这涉及6个步骤
1、提取组织或细胞总RNA后，除去占大部分的rRNA和tRNA，剩下编码RNA 和非编码RNA
2、对这些RNA进行测序，理想情况下，是直接检测，但是不现实，只有通过碱基互补配对的合成过程，才能知道原来样品中模板的序列，但是这个合成的长度是有限制的，所以只能先把这些RNA切割成小片段，再检测这些小片段的序列。具体原理见陈巍学基因视频。这个过程得到两种数据，一种是许许多多的碱基序列，一个是这些序列的表达频率。也就是一个是RNA是什么碱基序列，一个是RNA表达了多少量
3、由于上一步把RNA切割了，好像是一块拼图打散了，所以，这一步需要将这一个个的小块再重新拼成一个完整的图片。也就是比对，将检测到的RNA碱基序列，比对到参考基因组上，看某段RNA位于参考基因组的哪段序列上。这一步就好像一个拼好的拼图，上面有高高低低的小块，有些分子表达量高，它对应那个小块就高，反之就低。通过这一步，实验组和对照组都得到一个高高低低的拼图。
4、把实验组和对照组的拼图比较一下，看哪些RNA小块表达量是不一样的。或者你高我低，或者我高你低，从而得到这些差异表达基因名字的列表。因此这一步的结果都是一些基因名字或者转录本编号了。
5、将这些差异表达的分子，进行下游功能分析，比方看看它们都跟什么信号通路相关啊，可能跟什么功能有联系啊等等。这一步得到的就是很多结果图了。
6、下一步就是将筛选到的差异基因，结合你感兴趣的生物学功能或者过程，挑选出几个，进行再进一步的机制研究。这步就是湿实验了，也是决定文章层次的核心，这是需要人力和财力，再加上运气的事情，不过就算不做这一步，前5步也能发篇小文章灌灌水了。由于这一步涉及基础医学的机制研究方法，不在这篇文章的讨论范围内。
样品送测序仪器后，也就是上述第2步后得到会产生大量的数据，可能是多少个G的级别，有几万，甚至几十万的碱基序列，首先你要比对到参考基因组吧，然后你要看看实验组和对照组哪些基因表达有差异吧，其实这不是很复杂的事情，就是数据量太大了，如果就几十个，你完全可以用EXCEL查找，再标记，但是几十万个基因，谁能做到啊。所以现在需要一种工具，可以对数据进行批量编辑和操作。
感谢计算机发达的技术，前人早就帮我们想出来。Linux操作系统就可以实现对大量数据的批量编辑
，R语言可以实现大量数据的统计和做图。
好了，我们的下一步就是学习Linux操作系统和R语言了。
但是这两个部分包含了很多很多的知识，我们完全零基础，要是从头开始学，效率有点低，毕竟不是专业计算机出身，不需要一下子把所有东西都学会，先把目前需要掌握的学到，将来再举一反三，慢慢学其它的。
所以现在就开始模拟实战，从一个测序数据的样本开始，看看是怎么一步步得到文章中的结果的。
未完待续...
这两天宝宝得了幼儿急疹，耽误了几天，现在真的是上有老，下有小了，生活的压力会逼得你迅速成长起来。到这个时候才能深刻体会到时间是非常宝贵的，尽量少干不必要的事情，抓紧提升吧~
‘陆’ bulk RNA-Seq （4）合并表达矩阵
上一步定量，我们是对每一个样本进行的，现在需要将它合并到一个矩阵。这里通过一个脚本merge_to_matrix.pl 来合并。
输入：每个样本的定量结果
输出：
reads count 矩阵： genes.counts.matrix #原始矩阵
表达量矩阵：genes.TMM.EXPR.matrix #标准化后的矩阵
‘柒’ WGCNA（转载）
WGCNA原理及应用
WGCNA介绍：
WGCNA（weighted gene co-expression network analysis，权重基因共表达网络分析）是一种分析多个样本基因表达模式的分析方法，可将表达模式相似的基因进行聚类，并分析模块与特定性状或表型之间的关联关系，因此在疾病以及其他性状与基因关联分析等方面的研究中被广泛应用。
WGCNA算法是构建基因共表达网络的常用算法（详解： http://www.jianshu.com/p/94b11358b3f3 ）。WGCNA算法首先假定基因网络服从无尺度分布，并定义基因共表达相关矩阵、基因网络形成的邻接函数，然后计算不同节点的相异系数，并据此构建分层聚类树(hierarchical clustering tree)，该聚类树的不同分支代表不同的基因模块(mole)，模块内基因共表达程度高，而分属不同模块的基因共表达程度低。最后，探索模块与特定表型或疾病的关联关系，最终达到鉴定疾病治疗的靶点基因、基因网络的目的。在该方法中mole被定义为一组具有类似表达谱的基因，如果某些基因在一个生理过程或不同组织中总是具有相类似的表达变化，那么我们有理由认为这些基因在功能上是相关的，可以把他们定义为一个模块(mole)。这似乎有点类似于进行聚类分析所得到结果，但不同的是，WGCNA的聚类准则具有生物学意义，而非常规的聚类方法(如利用数据间的几何距离)，因此该方法所得出的结果具有更高的可信度。当基因mole被定义出来后，我们可以利用这些结果做很多进一步的工作，如关联性状，代谢通路建模，建立基因互作网络等。
WGCNA的用处：
这类处于调控网络中心的基因称为核心基因（hub gene），这类基因通常是转录因子等关键的调控因子，是值得我们优先深入分析和挖掘的对象。
在网络中，被调控线连接的基因，其表达模式是相似的。那么它们潜在有相似的功能。所以，在这个网络中，如果线条一端的基因功能是已知的，那么就可以预测线条另一端的功能未知的基因也有相似的功能。
下面的问答来自基迪奥，也能加深对WGCNA的理解
问1、调控网络和共表达网络有什么区别？
答：调控网络是个更广泛的概念，而共表达网络是调控网络的一种。
理论上我们可以利用各类信息构建调控网络（表达相关性，序列靶向关系、蛋白互作关系），另外调控网络构建的信息既可以来源真实的实验验证的关系，也可以来源生物信息的预测。而共表达网络特指利用基因间的表达相关性预测基因间调控关系的方法，而WGCNA又是共表达网络分析中最有效的方法之一。
问2、WGCNA分析适合的生物物种范围有规定么？
答：没有限制。对于任何物种中心法则都是存在的，调控关系对于任何物种都是存在的，所以WGCNA没有物种限定。
问3、同一物种，不同来源的转录组数据（比如不同文章/资料来源的），可以放在一起做WGCNA分析吗？
答：只要样本间有相似的生物学意义，是可以合并在一起做分析的。但要注意，不同批次之间的样本是有批次效应的，所以可能会带来一些误差，但是是可以放在一起分析的。
问4、相同材料不同处理之间，可以放在一起做WGCNA分析吗？比如重金属和盐碱处理。
答：可以的。这也正式WGCNA强大的地方，其可以将不同处理的样本，合并在一起做分析。其他方法则不一定有这么强大的能力，比如做基因表达趋势分析时，如果样本涉及到多个处理不同时期的时候，就不好合并分析（或合并后难以解读）。但WGCNA的方法关注的是调控关系，所以不管是多少个处理组，都可以很好的整合在一起做分析。
问5、不同批次的数据能放一起做WGCNA吗？
答：可以的。虽然有批次的干扰，但是干扰对WGCNA网络没有太大影响。因为WGCNA不是做差异分析，而是基因的共表达。因为批次效应理论上不影响相关性。
问6、不同类型的材料，比如亲本和F1，适合放一起进行WGCNA么？
答：如果是一个作图群体，当然亲本与F1是可以放在一起分析的，因为你只关心基因的表达模式，所以把亲本加进来是没有问题的。
问7、没有生物学重复，共3组，每组5个时间点能够做吗？
答：理论上有15个样本，是可以做WGCNA分析的。并且，分析出来的结果对你的研究应该是非常有用的。至少他会比趋势分析更有意义，更加准确。
问8、一般说WGCNA的样品不少于15个，15个样品考虑重复吗？不同倍性的材料呢？
答：15个样本这个是包含了生物学重复，比如5个时间点3个重复；在RNA-seq里面建议不要用不同倍性材料加进来。除非是有参考的多倍体，如果是无参的多倍体，不同倍性之间差异太大，会让调控网络不准确。所以用单一倍性的材料做调控网络会更加准确。
问9、可以将RNA-seq数据与蛋白组数据，甲基化数据放一起做WGCNA分析？
答：不能与蛋白数据一起分析。因为WGCNA是基于相关系数的算法。所以最好一起分析的数据变异度是类似的，RNAseq变异非常大，而蛋白的数据变异很小，两者的变化不在一个数量级上面。所以两种数据放在一起分析不合理。
但RNA数据可以尝试跟甲基化数据一起分析。当然我们也建议RNA数据与代谢组数据一起分析，因为代谢组的数据变异也非常大。
问10、表达量和表达的基因数目差异太大的样品可以一起分析吗？比如样品A有2k个gene表达 而样品B有2w个gene表达了 AB可以一起分析吗？
答：做WGCNA分析的时候，不能脱离生物学意义，既然要分析调控网络，那么应该分析有相似生物学意义的一组基因，比如说拿相似组织来一起做分析，比如不应该拿大脑的样本与脚趾的样本合并在一起做分析，因为很显然，这两个组织没有关联。如果两个样本之间是有相关联的生物学意义，哪怕表达的基因数不一样，或表达模式差异很大，那依然可以放在一起分析；但如果样本之间完全没有生物学意义，那么分析就没有意义。
问11、实验设计是case3个时间点（各点都有三个重复），control同样的3个时间点（每点三个重复），WGCNA怎么做？3个时间点和case-control两个因素能同时考虑进来分析吗？
答：可以的。做WGCNA是更加合理的，因为有两个梯度的样本，如果只是做差异分析的话，逻辑可能非常复杂，做WGCNA分析是对样本特性更好的解析，可以直观看到基因在六个处理组里面是怎样表达的。
问12、可以拿混合样本分析吗？比如一个病原细菌跟人类细胞的基因，能说明细菌跟人类细胞基因有调控关系吗？
答：可以。前提是病原菌有足够的数据并定量准确，并且这个分析是非常有意义的，最后可以说明这些病原菌可以调控哪些宿主基因。
问13、但是病原宿主混合分析的话，宿主蛋白不能分泌到宿主体内岂不是WGCNA生物学上也没有意义吗？
答：依然有意义。即使病原的基因没有分泌到宿主里面，但是病原的蛋白是会影响宿主基因的调控的，比如某个细菌感染某个植物，虽然细菌的蛋白不能直接分泌到植物体内，但会影响植物蛋白的分泌。混在一起分析依然是有意义，可以看到植物里面到底哪个基因对细菌蛋白产生应答作用。
问14、芯片数据两分类，每组20个样本，能否每组单独做WGCNA？
答：可以。WGCNA还有一种重要功能是做两个网络的比较，比如病人20个样本做一个调控网络，健康人做一个调控网络，然后两个网络做比较。
问15、WGCNA可以用来分析lncRNA对下游基因的调控分析吗？
答：可以。WGCNA网络有利于预测lncRNA的潜在功能。
问16、构建网络是用所有表达基因还是差异基因？
答：这个是具体问题具体分析。如果使用所有的基因分析，会导致运算量非常大。而也不是所有的基因在这个实验中都有生物学意义，所以我们会提前做一些过滤。
但用于分析的基因不一定是差异表达基因，有时可以用差异表达基因做一个并集，或通过计算变异系数将变异系数低的基因以及低表达的基因去除。但注意，如果你有关心的特定目标基因的话，应该尽量给予保留。
问17、关注某一个pathway上的基因以及调控因子之间的相关性，构建WGCNA网络的时候属于这个pathway的基因数量太少会不会影响结果呢？
答：这不是问题。在一个调控网络里面，样本的某个pathway上，并不是所有基因参与调控（或存在差异性），所以在做WGCNA分析的时候，会做一些过滤，将有变化的基因挑出来再做分析。即分析的是某个pathway上有变化的基因，不需要分析pathway上所有的基因，只需要分析那些变化的基因就够了。
问18、前期筛选的时候，要选出在所有样本中变异系数比较大的基因呢？还是直接用差异表达的基因取并集？用基因还是转录本，哪个好呢？
答：两则都可以，我推荐使用变异系数，选择那些变异较大的基因，来做下面的分析。然后建议用基因不要用转录本，因为转录本的定量是不准确的。
问19、变异系数一般取多大？
答：具体问题具体分析。例如，没有特定目标的时候，可以先计算变异系数，将变异系数的百分之前50来做分析，把变异系数偏低的后面一半过滤掉。
问20、输入数据用FPKM合适吗？
答：可以。
问21、RNA seq数据是RSEM值怎么办？
答：RSEM值原始输出结果为reads数，如果是RSEM值建议做一个RPKM校正再做分析。
问22、除了RPKM值以外，做WGANA是否还需要其他数据？TCGA数据可否来做WGCNA分析？
答：在做WGCNA分析必须要用表达量数据，但TCGA的数据某些层级没有表达量数据，没有表达量数据自然就无法做WGCNA分析。
问23、请问输入的基因样本的矩阵的时候，要不要对数据标准化？
答：做WGCNA分析的时候，不需要对数据进行标准化，输入RPKM值就足以做这个分析。虽然一些文章会做log2处理，但我认为取了LOG2后，会让一些表达关系没有那么丰富。
问24、每个样本有3个生物学重复，不需要对三个重复的表达量求平均值代表该样本吗？
答：注意，做WGCNA的时候每个样本是独立的，三个生物学重复样本是全部导入做分析，不是取均值再做分析，每个样本都是独立的。
问25、如果3个生物学重复，做WGCNA的时候是取三个值，还是用cuffdiff处理后取一个值？
答：如果是生物学重复样本进行调控网络分析，每个样本独立使用，而不是取均值。
问26、请问将样本信息同模块特征值进行相关性分析的时候，样本信息是怎么处理的呢？比如不同取样点、不同性别什么的，这不是数量性状信息的，这种情况应该怎么处理呢？
答：样本的任何信息都可以做模块相关性分析。比如相关时间点，可以按照先后量化为12134567。又如不同性别，男与女，可以定义为1，-1。任何性状量化为数字后，都可以进行相关性分析。
问27、怎么将模块与性状对应起来呢有些性状不好量化，如果直接将模块与分组对应，如何实现， 不需要量化指标么？
答：首先需要将性状量化，如果无法将性状量化，那么就无法分析。至于分组信息，也可以量化为类似00001111000（1代表一种组别，2代表另一组组别），实现分组信息的数字化。
问28、基因数量为3w左右时，moles数量为多少结果较为理想？怎么评价聚类效果的好坏？
答：moles数量没有标准，moles数量无法评估模块分的好坏，分组是否合理应该看树的树形图，比如树的分支很清晰就说明模块式清晰的。moles数量数由生物性状决定的。比如样本表达信息很丰富的时候，moles数量会很多；如果样本的基因表达相对单一，moles数量就会比较少。
问29、我运行例子的时候，得出来基因之间的direction全是undirected，这和前面的几种关系有什么区别？
答：WGCNA是一个undirected的方法，它的网络是无方向的，有相关关系但是无方向。
问30、如果做有向网络的构建，您推荐那些方法？
答：很多方法，例如贝叶斯的方法。
问31、非模式物种可以得出基因之间的相互关系类型么？得出的结果也是undirected么？
答：WGCNA是基于表达两处理的，所以即使是非模式生物，当然也可以他们之间关系，并且关系也是一个无向网络。
问32、选择几个表型数据进行结合分析比较好
答：越多越好，看实验设计。
问33、感染小鼠，5个时间点，3个重复，找不到合适的表型怎么办？
答：如果找不到合适表型，可以找某个时间点应答的基因，本身基因的表达趋势已经有某种生物学意义的。没有找到合适表型，也可以看变化趋势。不一定要做表型的相关分析，其他分析也是很有趣的。例如，可以对模块功能的富集分析，其实都是可以帮助你找到特定模块的。所以不用纠结于做某个表型的关联分析。
问34、weight就是tom值吗？
答：是的。
问35、剪模块是怎么做的？是根据TOM划分吗？需要自己设定，还是R自动的？
答：剪模块是R中自动完成的，不需要划分，但合并的时候你可以设定一个指标，比如差异度是0.25。
问36、看WGCNA说明是用相异矩阵D（D=1-TOM）去做聚类，然后动态剪切？
答：用TOM值来构建矩阵，TOM值就是两个样本的相似度，1-TOM值就是两个样本的差异度，相似度与差异度可以理解为一个东西，并不矛盾。
问37、模块特征值和样本性状相关分析的具体方法是？
答：R包用的是计算相关系数的方法。
问38、WGCNA里面一般会提到hubgene，如何确定hubgene?
答：在WGCNA分析里面，每个基因都会计算连通性，连通性高的就是hubgene。
问39、在R中安装“”WGCNA“”说不适合R3.3.1，那适合哪个版本？
答：WGCNA应该是所有版本都适合，如果版本没有可以考虑降低R软件的版本，这个对分析没有影响。因为不同R版本是一样的。
问40、用STEM分析的时候拟合多少个模型合适？
答：建议不要超过20个。模块太多不好分析。
参考网站：
http://tiramisutes.github.io/2016/09/14/WGCNA.html
http://www.jianshu.com/p/94b11358b3f3
http://www.omicshare.com/class/home/index/classdetail?id=20
‘捌’ 8.单细胞 RNA-seq 聚类分析：整合
目标：
挑战：
建议：
通常，在决定是否需要执行任何对齐之前，我们总是在 没有集成的情况下 查看我们的聚类。 不要总是因为您认为可能存在差异而总是执行集成 - 探索数据。 如果我们在 Seurat 对象中一起对两种条件进行归一化并可视化细胞之间的相似性，我们就会看到条件特定的聚类：
通常，在决定是否需要执行任何对齐之前，我们需要先查看没有整合的集群。 不要认为可能存在差异而执行整合——探索数据 。如果我们在Seurat对象中同时对这两个条件进行归一化，并可视化细胞之间的相似性，我们就会看到特定条件下的聚类:
不同条件下的细胞聚类，为了确保相同细胞类型的细胞聚集在一起，我们需要在不同条件下整合细胞。
为什么让相同细胞类型的细胞聚集在一起很重要？
在所有样品/条件/模式中我们要鉴别出 细胞类型 ，并在于我们的数据该细胞类型存在，我们需要探索我们的集群之间差异表达的基因。如果我们在一个簇中有我们的对照 T 细胞，而在另一个簇中有我们的受激 T 细胞，那么当我们试图在 T 细胞中找到富集的基因时，我们将返回许多作为 T 细胞标记的基因，以及相关基因处理和控制之间的差异。这些可能难以区分，使细胞类型识别更加困难。如果我们在同一簇中拥有相同细胞类型的所有细胞，它还可以通过识别 保守 的细胞类型来帮助我们识别细胞类型基因在簇之间差异表达。此外，如果您想对这些 T 细胞进行任何下游分析（例如控制/治疗之间的 DE 分析、配体-受体分析等），重要的是相同细胞类型的细胞存在于相同的簇中.
在本课中，我们将介绍不同条件整合样本，它改编自 Seurat v3 引导集成教程 。
如果细胞按样本、条件、批次、数据集、模式进行聚类，则此整合步骤可以大大地改进聚类和下游分析 。
为了整合，我们将使用来自每个组共有的高度可变基因（使用 SCTransform 识别），然后，我们将“整合”或“协调”这些组以覆盖具有相似“共同生物学特征”的细胞或组。例如，我们可以整合：
整合是一种强大的方法，它使用这些最大差异的共享资源来跨条件或数据集识别共享的亚群 [ Stuart 和 Bulter 等人。(2018) ]。整合的目标是确保一个条件/数据集的细胞类型与其他条件/数据集的细胞类型一致(例如，控制巨噬细胞与受刺激巨噬细胞一致)。。
具体地说，这种整合方法期望在组内的至少单个细胞子集之间有“对应”或 共享的生物状态 。集成分析的步骤如下图所示:
运行的不同分析步骤如下：
现在，使用我们的 SCTransform 对象作为输入，让我们执行跨条件的集成。
首先，我们需要指定我们要使用 SCTransform 识别的所有 3000 个最大突变的基因进行整合。默认情况下，此函数仅选择前 2000 个基因。
现在，我们需要 准备 用于整合的 SCTransform 对象 。
现在，我们将 执行 CCA，找到最好的伙伴或锚点并过滤不正确的锚点 。对于我们的数据集，这最多需要 15 分钟才能运行。 另外，请注意控制台中的进度条将保持在 0%，但要知道它实际上正在运行。
最后，我们可以 跨条件整合 。
整合后，为了可视化整合的数据，我们可以使用降维技术，例如 PCA 和 Uniform Manifold Approximation and Projection (UMAP)。虽然 PCA 将确定所有 PC，但我们一次只能绘制两个。相比之下，UMAP 将从任意数量的主要 PC 获取信息，以在这个多维空间中排列细胞。它将在多维空间中获取这些距离，并在二维空间中绘制它们，可以保持局部和全局结构。这样， 细胞之间的距离代表了表达的相似性 。如果您想更详细地探索 UMAP， 这篇文章 是对 UMAP 理论的一个很好的介绍。
要生成这些可视化结果，我们需要首先运行 PCA 和 UMAP 方法。让我们从 PCA 开始。
我们可以通过 PCA 映射看到 PCA 很好地覆盖了这两个条件。
现在，我们还可以使用 UMAP 进行可视化。运行该方法并绘制图表。
有时，如果我们 在 conditions 之间拆分绘图， 更容易查看所有细胞是否对齐，我们可以通过向 DimPlot() 函数添加 split.by 参数实现：
由于整合可能需要一段时间，因此 保存整合 seurat 对象 是一个很不错的选择。
‘玖’ RNA-Seq归一化问题
RNA-Seq在衡量基因表达水平时，若单纯以比对到基因上的reads数来计算表达量在统计学上是不合理的。影响因素有：
1.基因长度：需要基因长度来比较同一细胞内不同基因之间的表达。RNA-seq实验中众所周知的固有技术效果与基因长度有关：RNA（或cDNA）分子在测序之前先进行片段化，较长的转录本会比较短的转录本被剪切成更多的片段。因此，转录本的reads数不仅与其表达水平成正比，而且与其长度成正比。如下图，我们不能单纯的数比对到基因上的read数来比较表达量高低，需要考虑基因长度的影响。
这样来说，序列长的基因永远会被认为表达量较高，从而错误估计基因真正的表达量。为了消除基因长度产生的固有技术误差，在过去十年中，已针对RNA-seq数据开发了许多归一化方法，其中常用的有RPKM、TMM、RLE、upper quartile上四分位处理等。
所以需要对原始的表达矩阵进行标准化，去除掉测序深度和基因长度所带来的噪音。
是指所有细胞中均要稳定表达的一类基因，其产物是对维持细胞基本生命活动所必需的，管家基因是一类始终保持着低水平的甲基化并且一直处于活性转录状态的基因。
矫正的思路很简单，就是在变化的样本中寻找不变的量，那么在不同RNA-seq样本中，那些是不变的量呢？一个很容易想到的就是管家基因。但其实这种办法有一个非常强的先验假设：管家基因的表达量不怎么发生变化。其实管家基因有几千个，这几千个基因有一定程度上的变化是有可能的。因此这种方法不准确。
在RNA-Seq建库的过程中掺入一些预先知道序列信息以及序列绝对数量的内参。这样在进行RNA-Seq测序的时候就可以通过不同样本之间内参（spike-in）的量来做一条标准曲线，就可以非常准确地对不同样本之间的表达量进行矫正
数值概念：计算公式：CPM=C/N*1000000
设C为比对到geneA的read 数（read count ）,N为比对到所有基因的总read 数
用途：在某些情况下，只想了解每个基因被覆盖到的相对read数，而不希望对其做长度校正，就会使用这个指标。CPM 只对read count相对总read 数做了数量的均一化。如果想进行基因间表达量的比较，则不得不考虑基因长度的不同。如果进一步做长度的均一化，就得到了下面的RPKM.
RPKM，全称为reads per kilobase per million mapped reads，指的是每一百万个map 上的reads 中，map到外显子的每1K个碱基上的read数。
一般来说，基因越长，读取的次数（深度）越多，自然其有效读数就越多。而RPKM 就是为了消除这两个干扰的因素。以更好的比较不同结果。
数值概念：计算公式：RPKM=(1000000 C)/(N L/1000)
设C为比对到geneA的read 数（read count ）,N为比对到所有基因的总read 数，L为gene A的碱基数，RPKM法能消除基因长度和测序量差异对计算基因表达量的影响，计算得到的基因可直接用于比较不同样品间的基因表达差异
用途：用于与基因表达量相的后续分析，用于单端测序
计算步骤：
首先对总值数据进行标准化（当然正常来说map到的count 肯定不止这么多，这里只是除了10，但一般而言百万级的count 可以除一百万）
相比于RPKM，FPKM 计算的是fragments，也就是一对reads。与RPKM 的差别主要体现在，FPKM在一对reads map上的情况下只计数1，而RPKM 会计为2。适用于双端测序。
与 RPKM/FPKM 的差别在于，TPM 首先进行了基因长度的标准化，接着再进行了测序深度的标准化。
步骤：1.长度标准化
2.深度标准化
自然从这点来说TPM 的使用范围更为广泛。
参考 excel演示DESeq2归一化原理 - (jianshu.com)
要回答这个问题，我们需要先撇开所有形式上的计算，重新思考到底什么是RNA转录本的表达丰度这个问题。
事实上，对于任何一个制备好测序文库的待测样本，它上面任何一个基因的表达量（或者说丰度），都将已是一个客观存在的值，这个值是不管你改变了多少测序环境都不会变的。而且总共有多少个不同的基因在这这一刻进行了表达，实际上也已经是客观定好了
此刻，我们可以假定，对于样本X，其有一个基因g被转录了n次，同时样本X中所有基因的转录总次数假定是m次, 那么正确描述基因g转录丰度的值应该是：rg=n/m
同时，样本X中其他基因的转录丰度的计算也应该和上述公式类似。除了要把分子n换为其他基因对应的转录次数之外，分母m都一样。于是有趣的事情就是，所有基因转录本丰度的均值mean将是一个恒定不变的数，由以上定义这个数就是：mean = (a+b+c+...+n+...)/m/x
x代表基因个数
此时发现 (a+b+c+...+n+...)=m，公式化简变为mean = 1/x
这个期望值竟然和测序状态无关！仅仅由样本中基因的总数所决定的。也就是说，对于同一个物种，不管它的样本是哪种组织（正常的或病变的），也不管有多少个不同的样本，只要它们都拥有相同数量的基因，那么它们的mean都将是一致的。这是一个在进行比较分析的时候，非常有意义的恒等关系。
由于上面的结果是在理论情况下推导出来的，实际上我们无法直接计算这个r，那么我们可以尝试通过其他方法来近似估计r，只要这些近似统计量可以隐式地包含这一恒等关系即可
实际数据来证明FPKM和RPKM犯的错:
假定有两个来自同一个个体不同组织的样本X和Y，这个个体只有5个基因，分别为A、B、C、D和E它们的长度分别如下：
我们可以得到，样本X和Y的转录本的不变量，mean值是:
我们以FPKM的计算的为例子，以下这个表格列出的分别是样本X和Y在这5个基因分别比对上的fragment数和各自总的fragment数量：
于是，按照以上公式我们可以得到样本X和Y在这5个基因上的FPKM值分别为：
样本X在这5个基因上的FPKM均值FPKM_mean = 5,680;
样本Y在这5个基因上的FPKM均值FPKM_mean = 161,840
首先，我们可以把FPKM的计算式拆分成如下两部分。
第一部分的统计量是对一个基因转录本数量的一个等价描述（虽然严格来讲也没那么等价）:
第二部分的统计量是测序获得的总有效Fregment数量的百万分之一:
FPKM其实就只是这两部分的商！这有道理吗？分开来看它们似乎都有点道理，但是合起来的时候其实很没逻辑。
尤其是第二部分（N/10^6），本来式子的第一部分是为了描述一个基因的转录本数量，那么正常来讲，第二部分就应该是样本的转录本总数量（或至少是其总数量的等价描述）才能形成合理的比例关系，而且可以看出来FPKM/RPMK是有此意的，这本来就是这个统计量的目的。
可是，它却失败了！
N/10^6的大小其实是由RNA-seq测序深度所决定的，并且是一个和总转录本数量无直接线性关系的统计量——N与总转录本数量之间的关系还受转录本的长度分布所决定，而这个分布往往在不同样本中是有差异的！
比如，在有些基因中，虽然有效比对到基因的Fragment数是相等的，但一般来说长度越长的基因，其被转录的次数就越少。那也就是说，N必须将各个被转录的基因的长度考虑进去才能正确描述总体的转录本数！而FPKM/RPKM显然没有做到这一点，这便是FPKM和RPKM出错的内在原因。
‘拾’ RNA-seq比对到转录组
如果有注释较好的参考基因组的物种，可以直接将RNA-SEQ的数据比对到基因组上，也可以比对到从基因组注释出来的基因上，或者先比到基因组上再预测基因（也就是经典的tophat-cufflinks-cuffdiff流程）。
若没有注释较好的参考基因组，建议先用转录组序列从头组装（推荐trinity），再把从头组装的序列作为参考转录本，将RNA-SEQ数据比对到从头组装的序列中。或者把从头组装的序列和对应物种的NCBI unigene的序列整合在一起，去冗余之后再作为参考转录本。