1.Precision Medicine - exSEEK

Finding exRNA/cfRNA biomarker panel for cancer diagnosis
0) 背景介绍
我们都知道RNA是在细胞内转录产生的，由于RNA酶在环境中的广泛存在，在细胞外RNA很容易被降解。但是人们很早就意识到，在没有细胞的体液环境，如血清、唾液和尿液中，也可以稳定的检测到多种类型的 RNA(miRNA，mRNA, Y RNA,lncRNA,circRNA 等等)。这些存在于细胞外的RNA一般被称为 extracellular RNA(exRNA)或cell free RNA(cfRNA)。这些RNA通常不是以完全裸露的形式存在于体液环境中的，而是会由微囊泡 (microvesicles)、外泌体(exosome)包裹，或者与 RBP 密切结合形成 RNP 复合体。一般认为,这些载体(carrier)对exRNA/cfRNA有一定的保护作用，使得它们能在体液环境中免于降解。除此之外，一些结构稳定的RNA(or a RNA domain)会因为自身二级和三级结构可以抵抗 RNAase 的降解， circular RNA因为没有游离的5'和3'末端也会对外切酶的降解有比较强的抵抗作用。
exRNA/cfRNA起初是来自于细胞内的，所以exRNA的表达谱和转录后修饰等在一个侧面上反映了人的生理状态。血清、唾液和尿液都是容易通过无创的方式获取的样本，所以体液exRNA/cfRNA的一个重要应用场景是作为一类无创检测的标志物，服务于人体健康状况检测和疾病的诊断，如孕检，癌症的早期诊断，肿瘤的监测和预后等。
目前最常见的对exRNA/cfRNA进行高通量检测的实验手段就是我们前面反复介绍的二代测序。通过前面的学习，我们已经对二代测序数据的处理有了一定的了解。本 Quiz希望读者利用生物信息学工具从二代测序数据中提取特征，使用一些机器学习方法，评估exRNA/cfRNA对于健康人和癌症患者的区分能力，并给出一个有预测能力的特征组合。
1) 数据集介绍
我们提供了血浆long RNA双端测序和small RNA单端测序产生的两个数据集，以下是数据分析的整体流程:
我们已经对原始数据进行了预处理(去除adapter、去除低质量序列等)和mapping，请大家直接从两个数据集的bam文件(bam文件已按坐标sort过)出发进行后续分析。
相关的数据都放在P集群的/data/2022-bioinfo-shared/data/quiz-I目录下。
数据集1 (long RNA)
该数据集包含373个样本，样本来源于结肠癌(CRC)，胃癌(STAD)，肺癌(LUAD)，食管癌(ESCA)和肝癌(HCC)5种癌症患者，以及健康人(HD)。
数据位于/data/2022-bioinfo-shared/data/quiz-I/exRNA-long目录。文件内容解释如下:
├── bam
│   ├── CRC-2124325.bam
│   ├── CRC-2124325.bam.bai
│   ├── CRC-2211756.bam
│   ├── CRC-2211756.bam.bai
...
├── genome
│   ├── bed
│   │   ├── lncRNA.exon.bed
│   │   ├── lncRNA.intron.bed
│   │   ├── mRNA.exon.bed
│   │   ├── mRNA.intron.bed
...
│   └── gff 
│       └── gencode.v38.annotation.gff3
├── metadata.txt 
└── README.md
bam目录下是每个样本的reads向基因组mapping得到的bam文件(按坐标排序)和bam index，每个样本对应一个bam文件
genome目录下是需要用到的基因组注释信息
metadata.txt中提供了样本的相关信息，示例如下:
sample_id	label	source	sex	library	dataset
CRC-2124325	CRC	P1	M	reverse	train
CRC-2211756	CRC	P1	M	reverse	train
CRC-2276341	CRC	P1	F	reverse	test
CRC-2277930	CRC	P1	M	reverse	train
metadata一共包含6列，每列的含义解释如下
sample_id: 样本id，每个样本的id是唯一的
label: 样本的标签，也是我们实际预测的对象
source: 样本来源的医院
sex: 样本供体的性别
library: 建库的链特异性，有reverse/forward两种
dataset: 参考数据集的一种人为划分，我们把一部分人的数据用来做training， 另一部分人的数据用来做 test。注意：对于本测试中的这种较小的数据集，我们不推荐用这种人为一次性划分 traning 和 test 的方式。更好的方法参考下文的 ”数据集划分“。
数据集2 (small RNA)
公共数据GSE71008​
包括192个样本，来源于结肠癌(CRC)，前列腺癌(PC)和胰腺癌(PAAD)以及HD。
数据位于/data/2022-bioinfo-shared/data/quiz-I/exRNA-short目录。文件内容解释如下:
.
├── bam
│   ├── Sample_1S1
│   │   ├── lncRNA.bam
│   │   ├── lncRNA.bam.bai
│   │   ├── miRNA.bam
│   │   ├── miRNA.bam.bai
...
│   ├── Sample_1S10
│   │   ├── lncRNA.bam
│   │   ├── lncRNA.bam.bai
│   │   ├── miRNA.bam
│   │   ├── miRNA.bam.bai
...
├── genome
│   ├── bowtie2-index
│   │   ├── lncRNA.1.bt2
│   │   ├── lncRNA.2.bt2
...
│   └── fasta
│       ├── lncRNA.fa
│       ├── miRNA.fa
...
└── metadata.txt
注意 small RNA 数据和上述的 long RNA 数据的处理方法不同，bam目录下是每个样本的reads向转录组mapping得到的bam文件(按坐标排序)和bam index，所以对每个样本都会对应9个bam文件(miRNA,lncRNA,mRNA,piRNA,snoRNA,snRNA,srpRNA,tRNA和YRNA各自对应一个bam文件)。详见下文的“Different Mapping stratgies for small and long RNA-seq data”
genome目录下是我们在mapping时使用的序列和对应的bowtie2 index
metadata.txt中提供了样本的相关信息，示例如下:
sample id	RNA Isolation batch	library prepration day	gel cut size selection	label	sex	id	stage	dataset
Sample_1S1	2	22	7	CRC	F	GSM1825212	I	train
Sample_1S10	3	24	1	CRC	F	GSM1825215	I	train
Sample_1S11	3	24	2	CRC	F	GSM1825239	I	test
Sample_1S12	3	25	7	CRC	M	GSM1825071	I	train
metadata一共包含9列，每列的含义解释如下:
sample id: 样本编号，对于每个样本是唯一的
RNA Isolation batch, library prepration day, gel cut size selection三列记录了实验批次的信息
label: 样本的标签，也是我们实际预测的对象
sex: 样本供体的性别
stage: 肿瘤分期的信息
id: GEO数据库中的样本id
dataset: 参考数据集的一种人为划分，我们把一部分人的数据用来做training， 另一部分人的数据用来做 test。注意：对于本测试中的这种较小的数据集，我们不推荐用这种人为一次性划分 traning 和 test 的方式。更好的方法参考下文的 ”数据集划分“。
Different mapping stratgies for small and long RNA-seq data
long RNA-seq: 长RNA双端测序的mapping和RNA-seq的常规分析流程比较一致，即用spliced aligner(我们这里使用的是STAR)将reads比对到参考基因组hg38上。所以对于长RNA数据,一个样本对应一个bam文件,bam文件的record都在hg38的坐标上。
small RNA-seq: 小RNA数据集的reads长度比较短，如果直接往基因组上比对容易产生大量的multiple mapping。所以我们这里采用的方法是人为的对不同类型的RNA定义一个优先级，按这个优先级先把reads往前一类RNA的转录本序列上mapping(DNA-seq aligner即可，这里用的是bowtie2)，再把unmapped reads往下一类RNA的转录本序列上mapping。这种"sequential mapping"是处理小RNA测序数据的一种常见策略，大家可以参考2019, Cell Systems, Joel Rozowsky et al., exceRpt: A Comprehensive Analytic Platform for Extracellular RNA Profiling这篇文章。优先次序的定义有一定的主观性，通常会给注释的可信度比较高的RNA类型和我们在研究中比较关注的RNA类型较高的优先度，我们这里采用的顺序是miRNA,lncRNA,mRNA,piRNA,snoRNA,snRNA,srpRNA,tRNA,YRNA。对于每一个小RNA测序样本，我们都会提供多个bam文件，对应不同的RNA类型，每一个bam文件的record都在相应的转录本序列的坐标上。
2) 报告要求
提交一份完整的工作报告(中英文不限)，提供主要的代码。
请读者使用我们提供的数据，参考 Help 中提供的相关信息，完成以下工作:
Part I. Prepare count matrix
RNA-seq最直接的作用是用于基因表达的定量，利用RNA测序数据估计出有注释的基因或转录本的相对丰度是一个很容易想到的，可能在癌症病人和健康人之间有区分能力的特征。在part I中，我们要求大家对两套RNA-seq数据，分别构建表达矩阵。
long RNA:
用featureCount一类的工具统计每个基因的read counts。
根据bam文件计算在我们考虑的RNA中，不同RNA类型(例如可以按mRNA,lncRNA,snoRNA,snRNA,srpRNA,tRNA,YRNA,pseudogene和intron来划分，请参考help中的相关介绍)所占的比例。
把不同样本的count合并到一个count matrix中。
注意：
我们这里用到的长RNA数据都是链特异性建库产生的，但是这些数据当时在做实验时分别来自2个建库方式：有一些是forward stranded，另一些则是reverse stranded。这一信息在metadata.txt中已标明，在count时需要使用相应的参数。
small RNA:
我们提供了一个脚本scripts/count-transcripts.py用来统计每一个转录本上有多少reads。在基本的分析中我们建议大家只考虑miRNA和piRNA，更多的分析可以参看下面的 "Part IV Optional Tasks".
把不同样本的count合并到一个count matrix中。
Part II. Matrix processing
在该部分中，我们要求大家对count matrix进行预处理，将矩阵处理成机器学习模型直接可以接受的形式。我们要求大家对矩阵进行如下处理：
1.过滤低表达的基因
2.对表达矩阵进行normalization
3.用计算方法去除batch effect
4.通过数据可视化，分析预处理有没有达到预期的效果
5.对特征进行scaling(如果只使用树模型该步骤可省略)
Part III. Machine Learning
long RNA
基本分析: 癌症患者(不区分癌症类型) vs. 健康人(healthy donor, HD)的二分类问题
进阶分析: 尝试特定一种癌症类型和健康人的二分类问题，以及不同癌症类型间的多分类问题
small RNA
基本分析：CRC vs. HD之间的二分类问题
进阶分析：在男性样本中尝试PC vs. HD，以及CRC vs. PC的二分类问题(女性不会患前列腺癌)
数据集划分
我们用到的两个数据集都比较小，只做一次训练集和测试集的划分，在测试集上的表现可能会有比较大的随机性，更合适的做法是通过多次random split, 或多次k fold cross validation，或bootstrap的方式产生多个训练集和测试集，再用每一对训练集和测试集重复同样的计算(在训练集上进行调参，特征选择和模型拟合，在测试集上进行性能评估)。所以我们提供的数据计划分仅作为参考，有精力的同学可以尝试多次划分训练集和测试集，分析模型表现的分布。
特征选择
请大家用任意一种方式对相应的分类问题进行特征选择，鼓励大家在实现特征选择的过程中考虑更多的实际中存在的问题，或使用更复杂的特征选择策略。请大家根据自己的时间和兴趣选择和尝试下面的不同策略：
很多模型自带特征选择的功能。所以最简单的做法是直接用这样的模型(如LASSO)在训练集上拟合模型，交叉验证选择超参数，汇报系数不为0的特征。在1.3.feature-engineering一节中我们把这类方法称为"embedded"方法。这样得到的模型可以直接用来在测试集上进行性能评估。除了基于"embedded"方法的特征选择，我们还希望大家多做其他尝试，如前面介绍的"filter"和"wrapper"等方法。
特征选择时还可以引入先验知识。例如2021, Nature Communications, Matthew H. Larson et al.专门挑选一些在所有的健康人中都几乎检测不到的基因，文中称为_dark chanel biomarkers_，希望以此提高诊断的特异性。此外，有的研究也将特征选择的范围限制在一些已知的“Cancer Genes”上。
对冗余特征的处理: 生物学上有些基因是共表达的，计算上也有可能引入一些人为的相关性。大家可以考虑在特征选择前就去除一些高度相关的特征，或筛选出特征后去除高度相关的特征。
评估不同特征选择方法的稳健性: 从训练集sample一部分样本(例如80%的样本)用于特征选择，重复多次，分析每次选到的特征的重合情况。有很多指标可以用来描述特征选择的稳定性，如Kuncheva index等。
注意:
特征选择原则上和模型拟合一样，只能用 training set 的信息，不能用到test set的label。
模型评估与特征解释
1.用热图对选中的特征在训练集和测试集上的数值大小进行可视化
2.绘制测试集上的ROC曲线，计算AUROC
3.汇报挑选不同数量的feature时分类效果，用AUROC作为指标绘制折线图
4.分析挑选出的feature的生物学意义
该项目的重点和目标在于特征提取与特征选择，最后形成一个希望能用于诊断的RNA特征组合。可靠的，可解释的特征比测试集上准确性的数值大小更为重要。所以建议大家在特征的选择和解释上多花一些精力。
Part IV. Optional Tasks
基于RNA-seq数据，除了 count gene expression 构建 expression matrix 这种常规方法。我们还可以做很多其他分析和探索，大家可以计算一些新的feature重复Part II和Part III，并分析这些特征的分类效果。为了帮助大家开拓思路，我们给出几个可能的方向：
Small RNA: peak calling
对小RNA数据我们可以不局限于常见的 miRNA 和 piRNA， 虽然小RNA测序数据中测到的reads主要来自miRNA，piRNA等小RNA，但还有些其他的：
1.Peak Calling inside the long RNAs: 有些 small RNA reads来源于长RNA（如mRNA，lncRNA）的片段（long RNA derived fragments）。这些fragments可能多数来源于一些长RNA在血浆中比较稳定的，不容易被降解的局部区域。我们可以用计算方法(如peak calling等)定义出这些区域，再选择性统计这些区域(而非全长转录本)的reads数量(在前面的流程图中我们将其称为"domain features")。这种做法一方面在生物学上相对更加合理，另一方面也有希望减少噪声的影响，提高分类效果。在 Help 部分，我们会对具体做法给出更详细的说明。已有研究者对此进行了一些研究，请大家参考：
2020 elife - Identification of protein-protected mRNA fragments and structured excised intron RNAs in human plasma by TGIRT-seq peak calling (2020, eLife, Jun Yao et al.)
2.Peak Calling in the unannotated regions: 还有些 reads来源于一些过去没有注释的区域，比如之前的一篇研究发些了一些新的 RNA，称为 oncRNA (orphan noncoding RNA)。相关文章：
2018 Nature Medicine - Cancer cells exploit an orphan RNA to drive metastatic progression
2021 Gut - Unannotated small RNA clusters associated with circulating extracellular vesicles detect early stage liver cancer
Long RNA: explore more genomic regions
对长RNA数据我们可以不局限于一些常见的RNA的exon，我们还可以考虑一些在gencode v38中没有被注释成常规exon的区域，如一些重复序列(TE: Transposable Elements)和其他未注释的转录本（如上面的 oncRNA）等，还可以考虑 circRNA 等，分析align到这些区域的reads数量对于癌症患者和健康人有没有区分能力。更多的其他 Genome Regions可以参考 Appendix VI. Genome Annotations.
Other variances derived from RNA-seq
我们还可以从RNA-seq数据中计算RNA splicing, RNA editing等其他变异。
休息一会
Grail
随着科学技术的不断发展，尤其是 21 世纪初高通量测序技术（NGS） 的出现，使肿瘤诊断从传统的病理和影像学检测跨入精准诊断时代，“液体活检”的概念也应运而生。液体活检（Liquid Biospy）是一种利用 Sanger、qPCR、NGS 等基因测序技术从血液、脑脊液、唾液等非实性生物标本中检测循环肿瘤 DNA（ ctDNA） 、循环肿瘤细胞（ CTCs） 、外泌体（exosomes） 等生物标志物的肿瘤诊断方法。
液体活检作为可用于癌症早筛的一种无创检测技术，一直以来备受科研和临床研究的关注。测序巨头Illumina首席执行官Jay Flatley此前在接受媒体采访时曾表示，“液态活检”的市场规模至少达400亿美元，甚至宣称这项技术可能是癌症诊断领域最激动人心的突破。
GRAIL正是一家以“液体活检”为中心的公司，被外界称为全球癌症血液筛查公司中的领先者、癌症大数据领域的独角兽 。2016年1月由基因测序巨头Illumina联合比尔盖茨、Bezos Expeditions、和Sutter Hill Ventures等投入一亿美金，并且分出一部分公司骨干成立了Grail。据统计，GRAIL自2016年成立以来，短短两三年时间就已经获得了全球多家公司超过15亿美元的巨额投资，它成为历史上融资规模最大的三家生物技术公司之一，包括腾讯在内的一些中国公司也是其资方。
Jennifer Cook是GRAIL的目前的首席执行官。此前，Jennifer曾在Roche Pharmaceuticals / Genentech担任过多个高级管理职位，负责产品开发和商业化的整个生命周期。2016年 ，因其对医疗保健行业的贡献和鼓舞人心的领导力而获得认可，她被医疗保健女企业家协会评为年度女性。Jennifer拥有斯坦福大学人类生物学和生物学硕士学位，以及加州大学伯克利分校哈斯商学院的MBA学位。
最近，在2018年欧洲肿瘤内科学会（ESMO）年会上，GRAIL公司发布了CCGA（ Circulating Cell-free Genome Atlas ）研究项目的最新数据。当前研究结果显示，利用血液进行癌症早期筛查不仅可行，而且在不同类型癌症中还具有高度特异性。使用无创、简单、精准的液体活检方法进行癌症早期筛查，代表着人类征服癌症的最大希望。现在，这个梦想变得越来越可触可及了。
Last modified 9mo ago