人类结直肠癌单细胞多组学分析
文献详解栏目
每个人的时间精力有限，必须优先阅读相关文献，开设这个栏目也是希望为大家推荐高质量的单细胞相关文献。如果大家对单细胞转录组感兴趣可以关注一下，哪怕每天只学一点点，积土成山，积水成渊。
当然一个人的力量终归是小的，我也希望汇聚一群人，形成一个场，这里头最重要的生产力不是单个人多聪明，多厉害，而是每个人相互作用，形成的那个氛围。
希望大家能有所收获！
文章信息
文章由北医三院付卫团队、乔杰团队和北大生命科学学院汤富酬团队于2018年11月合作发表在Science上，文章题目是：Single-cell multiomics sequencing and analyses of human colorectal cancer
这一篇理解难度稍大，当做背景知识的了解，至于其中数据的分析细节文中没有提及
导语
之前一次只能研究单细胞层面的基因组、转录组或DNA甲基化组其中一种，不能在一个细胞中同时研究多个组学，2016年汤富酬研究组将三重组学研究方法scTrio-seq（single-cell triple omics sequencing technique）发表在Cell Research上。2018年11月北医三院付卫团队、乔杰团队和北大生命科学学院汤富酬团队又在Science合作发表了Single-cell multiomics sequencing and analyses of human colorectal cancer
文章优化了单细胞多组学测序技术，并对原发瘤、淋巴转移和远端转移区域分别采样，首次从单细胞分辨率解释了结直肠癌的发生与转移过程中中基因组拷贝数变异、DNA甲基化异常及基因表达改变的特点，证明了用单细胞多组测序重建遗传谱系和追踪其表观基因组和基因表达动力学的可行性
背景知识
淋巴转移（lymphatic metastases）和远端转移(distant metastases)
癌症可以以2种方式出现在淋巴结中：一种是从淋巴结形成的肿瘤叫淋巴瘤，另一种是从其他部位扩散叫淋巴转移（更为常见）。淋巴转移一般会和乳腺癌、前列腺癌、肺癌、结直肠癌的不良预后相关，淋巴转移虽然不是致死因素，但会导致癌细胞扩散到重要器官。远端转移也叫恶性转移，肿瘤细胞从原始发生的部位借由侵入循环系统，转移到身体其他部位继续生长，几乎不可能使用外科手术切除根治，多半只能用大范围循环全身的放射治疗或化疗等手段来抑制已转移的癌细胞。
2017年发表在Science的Origins of lymphatic and distant metastases in human colorectal cancer 中描述了结直肠癌之前的肿瘤扩散的TNM层级是：primary tumor(T)=》lymph node system(N)=》distant metastases(M)，但是有临床证据表明移除淋巴结并不会提高病人存活率，因此N和M之间的关系可能并不是简单的上下级。文章发现淋巴结和远端器官存在独立起源的证据，只是淋巴结转移机制形成的更快，因此淋巴转移形成更早，发生更频繁。因此作者不推荐直接假设淋巴结会引发远端转移而直接切除(https://www.medscape.com/viewarticle/882502)
结直肠癌colorectal cancer(CRC) ：是结肠癌和直肠癌的统称，是消化道恶性肿瘤之一。2018年Cancer Stats统计显示：结直肠癌在男性中发病率第2死亡率第3，女性发病率第4，死亡率第3。
约95％的结直肠癌是由结肠和直肠内壁的腺细胞发展而来，癌症通常开始于内壁最内层，并缓慢生长到外层。(http://www.cancer.org/Cancer/ColonandRectumCancer/DetailedGuide/colorectal-cancer-what-is-colorectal-cancer)
结直肠癌发生率在40岁开始增加，60～75岁达高峰，结肠癌在女性患者较常见；直肠癌在男性患者常见；大约5%的结肠癌或直肠癌患者在结肠和直肠有两个或更多病灶，并非简单从一个病灶转移至另一个所致(https://www.msdmanuals.com/)
结直肠癌典型的分子特征是：基因组不稳定性、表观遗传学异常、基因表达紊乱
瘤内异质性 Intratumoral heterogeneity (ITH) ：恶性肿瘤的特征之一，肿瘤异质性包括肿瘤间异质性（不同肿瘤细胞之间的基因与表型不同）和肿瘤内异质性（相同肿瘤细胞以内的基因与表型也不同），其中肿瘤内异质性又包括空间异质性（相同肿瘤不同区域不同，如未扩增的细胞背景中有成簇扩增细胞；少量扩增背景中有未扩增的细胞；孤立的细胞扩增【利用多位点取样方案或者tissue microarrays (TMAs) 调查】）与时间异质性（原初肿瘤与次生肿瘤不同）。
异质性的产生是因为同一肿瘤由多种不同基因组特征的细胞组成，每一种细胞构成一个亚克隆(subclone)。肿瘤组织会存在对治疗药物有抗性的亚克隆，但比例不高。当治疗的药物除去敏感的亚克隆时，抗性的亚克隆细胞会不受药物抑制并且少了空间竞争，因此会加快生长速度，导致肿瘤复发或者发生转移，而且转移后的亚克隆对同种治疗方案也会有抗性。因此研究肿瘤细胞的亚克隆以及不同的亚克隆的转移是一个热点，尤其是亚克隆是如何从原位癌转移到其他脏器而形成转移癌。2017的一篇文章专门研究了转移癌亚克隆与原位癌亚克隆的进化关系 ，他们发现结直肠癌肿瘤转移癌高深度测序就可以找到肿瘤的大部分基因组变异，另外验证了"转移癌多克隆起源说(转移癌是由多个起源于原位癌的亚克隆发展而来，而非由单个细胞发育)"，发现了结直肠癌的淋巴与远端并行转移 。
方法
scTrio-seq2技术：用于somatic copy number alterations (SCNAs)拷贝数变异、DNA甲基化特征、细胞连续的转录信息；整合了单细胞重亚硫酸盐测序（scBS-seq）用于全基因组甲基化分析；研究的细胞数量从之前的25个增至1900个
分析了12个CRC患者（III期或IV期）的约1900个单细胞，7.6Tb高质量测序数据。DNA甲基化研究中平均每个细胞测序量为4.1Gb，平均覆盖到全基因组内870多万CpG位点；转录组研究中每个细胞测序量为0.9Gb，平均覆盖3700多个基因
多区域采集了10个患者的原发瘤、淋巴转移瘤或远端转移瘤样本(利用两种不同来源的细胞，可以发现每个患者因突变而产生的遗传谱系)
文章实验流程
图A是取样：治疗前后的肿瘤区域（包括原发瘤primary tumor，PT；淋巴结转移位 lymph node metastasis, LN；肝转移位 liver metastasis, ML；化疗后肝转移位 posttreatment liver metastasis, MP），然后测序分析了基因组、转录组、甲基化组；
图B是化疗6个周期后的患者CRC01取样：一共取了ML（4个）、MP（5个）、LN（3个）、PT（4个）共16个肿瘤区域
单细胞甲基化数据处理
首先raw reads去接头、引物、低质量碱基，然后利用
Bismark（V0.7.6）clean reads比对到hg19基因组，PCR重复利用samtools rmdup(V 0.1.18)去除，数据统计(比对数、比对率、CpG位点数、亚硫酸氢盐转化率等)[其中亚硫酸氢盐转化率是由非甲基化的lamda DNA的spike-in计算的]，CpG位点小于200万个或亚硫酸氢盐转化率小于98.5%的细胞被排除。仅使用甲基化水平大于0.9或小于0.1的CpG位点进行总体甲基化水平计算。启动子区设定为转录起始位点的上游1 kb到下游0.5 kb。为了计算RefSeq基因各基因体的DNA甲基化水平，将每个基因体划分为100个等分，并将其上下游侧翼区域(15 kb)分别划分为10个等分。基因组注释信息从UCSC获取。利用bedtools(V 2.17.0)和自定义脚本(没放Git链接)计算平均甲基化水平，设置滑动窗口大于等于3个CpG位点
根据甲基化测序数据估计拷贝变异数
主要基于Garvin等人开发的Ginkgo算法 ，基因组被分成10856个不等长的bins，长度中位数为250kb，并根据算法的过滤器排除了一些异常的bins。
BED文件是利用bedtools从BAM文件得到。每个bin的read counts值利用所有bins的count平均值进行标准化，采用低水平均一化(Lowess normalization)来校正基因组GC含量的偏差。另外，以正常的二倍体细胞作为对照，减少scRS-seq的其他误差。利用Circular binary segmentation (CBS) 对copy number文件进行分隔，参数为"
alpha = 0.0001"和"undo.prune = 0.05" 。CBS分隔后，每一段的所有bins的计数重置为这一段的bin count的中位数。每个单细胞的基本拷贝数由smallest sum-of- squares (SoS) error和 scaled copy-number profile (SCNP)决定，其中SCNP又进一步四舍五入，最后得到了整数值copy-number profile (FCNP)。利用 GISTIC2.0 (https://software.broadinstitute.org/cancer/cga/gistic) 鉴定了重要的SCNAs和潜在的基因靶点。TCGA数据分析
从https://tcga-data.nci.nih.gov/docs/publications/coadread_2012/获得已发表的人类CRC的SCNA片段数据，用于SCNA频率统计(不包括X染色体)。将本文的数据与TCGA的进行比较时，本文研究的CMS3类型的患者CRC02是找不到对应的，只有Affymetrix SNP 6.0芯片得到的178个non-hypermutated样本。SCNA谱进一步转变成长度不等的bins用于单细胞SCNA统计，拷贝数为>2.5的bin表示扩增，小于1.5表示缺失。Circos图是根据https://github.com/venyao/shinyCircos制作的。
WGS数据处理
raw reads =》trimmed =》BWA mem(V 0.7.12) 比对到hg19 =》 samtools sort =》Picard(V1.139) merge BAM文件 + 标记重复 =》BAM文件用GATK(V 3.4-46) 预处理=》muTect (V 1.1.4) call SNVs ，自己脚本过滤=》取每个患者的外周血或邻近正常组织作为对照(somatic variants)，在线Venn图做出SNV数量
单细胞RNA-seq处理
利用汤教授自己的方法得到的数据处理：预处理过程都一样，然后用STAR(V 2.5.0) 2步比对到hg19，Cufflinks(V 2.2.1) 使用默认参数进行FPKM定量；
利用multiplexed scRNA-seq方法得到的数据，先利用read2的barcode信息将reads分配到每一个细胞，每个细胞中read2对应的read1利用read ID分隔，read1中的TSO序列 利用自己的脚本过滤掉，然后利用Tophat(V2.1.1) 单端比对到hg19，利用UMI实现TPM标准化(文中说道：大部分的表达量都使用
log2(FPKM + 1)或者log2(TPM/10 + 1)) ，然后统计了mapped read numbers, mapping ratios, RefSeq gene numbers等信息，根据FPKM > 1 or TPM > 1去除了比对率 < 20%或有表达量的Refseq基因数量 < 1500
结果
根据甲基化数据得到的基因组拷贝数变异来推断基因谱系
Genomic alterations in tumors provide markers for lineage tracing. 克隆变异出现在肿瘤早期阶段，亚克隆拷贝数变异标志着亚型的出现
结直肠癌患者单个癌细胞的染色体拷贝数变异谱+高精度的染色体内断点信息=》谱系追踪=》原发位肿瘤（PT）的亚克隆结构通常比转移位肿瘤更复杂
结果得到了5个患者的90个细胞以上的甲基化数据，细胞被分成了不同的基因亚型
其中，CRC01基于21个亚克隆的拷贝数断点，鉴定了来自2个不同谱系（A、B）的12个亚型，其中每个亚型都有4-8个亚克隆的断点(断点的上下位置和拷贝数变异数增加、减少对应)，A5亚型同时出现在了肝转移位和淋巴转移位，表明这两种移位有共同起源，这5个病人的原癌亚克隆结构比其他转移类型更复杂
CRC01患者癌细胞的单细胞染色体拷贝数变异谱：
CRC01患者亚克隆拷贝数断点：
CRC01患者亚克隆结构：
甲基化异质性
结肠直癌细胞的DNA甲基化水平要低于癌旁的正常上皮细胞，同一肿瘤组织中同一谱系的甲基化程度相近，不同谱系出现差别。低甲基化基因组区域显著富集在LTR（long terminal repeats）、LINE-1 (long interspersed nuclear elemnt 1)和异染色质区域(H3K9me3)，而高甲基化的基因组区域显著富集在CpG岛、H3K4me3和开放染色质区域。
以CRC01患者为例：
甲基化的异质性主要来自同一个患者肿瘤内不同亚克隆之间的DNA甲基化差异，而不是同一个亚克隆内部不同细胞间的DNA甲基化差异
甲基化谱和基因表达谱的相互关系
启动子区域的甲基化与相应基因的表达呈显著的负相关，而基因区的甲基化与相应基因的表达呈正相关
【The gray lines represents individual cells. The blue line represents the mean value for each patient. TSS, transcription start site; TES, transcription end site.】
同一遗传谱系的肿瘤细胞在转移过程中DNA甲基化及基因表达的变化情况
同一患者同一谱系的肿瘤细胞从原发灶到转移灶全基因组DNA甲基化水平基本稳定，组内局部区域可能会有比较大的波动
结直肠癌去甲基化特点
每个亚型内的去甲基化 程度是一致的，不同亚型程度不同 ；
正常上皮细胞的基因组区域甲基化越高，它就越容易发生去甲基化；
去甲基化的程度与基因组重复序列L1（long interspersed nuclear elemnt 1）以及癌旁正常组织中H3K9me3修饰的密度呈正相关，与H3K4me3标记和正常组织的开放染色质区域密度呈负相关；
有趣的是，L1比LINE-2更活跃，在所有病人的癌细胞中显示了更强的去甲基化能力，这个与胚胎发育中情况相反(胚胎发育过程中L1一般比L2去甲基化能力弱) ，说明在肿瘤发生与发展过程中，L1和异染色质区域产生了异常的去甲基化过程，打破了正常的发育规律
癌细胞相比于癌旁细胞的DNA甲基化水平：
染色体水平的DNA甲基化异常和基因组不稳定性
结直肠癌细胞中6条染色体（4号、5号、8号、13号、18号、和X染色体）倾向于发生更强烈的DNA去甲基化，其中三条低甲基化染色体（8、13和18）在TGCA和研究的患者中都有较高的拷贝数变异。结合WGS结果发现，有5条第甲基化的染色体(4号、5号、8号、13号、和X染色体)的单核苷酸变异(SNVs)发生显著富集