1 各种定义

1.1 驱动基因

与癌症发生发展相关的重要基因称为驱动基因，驱动基因决定了这个癌症的最主要的原因。当驱动基因突变后，就会把癌细胞“驱动”起来。

1.2 mutation rate of each sample

calculated by number of mutations per megabases (Mb) of target sequences.

1.3 clonal hypothesis

suggested that multifocal tumors arose from the same origin, a founder tumor, and the other tumors were formed by intraluminal or intraepithelial seeding of cells from the founder tumor.

1.4 field hypothesis

emphasized the significance of tumor microenvironment in multiple tumorigenesis.

1.5 SNP

据估计，人类基因组中每1000个核苷酸就有一个SNP，人类30亿碱基中共有300万以上的SNPs。SNP 遍布于整个人类基因组中，根据SNP在基因中的位置，可分为基因编码区SNPs（Coding-region SNPs，cSNPs）、基因周边SNPs（Perigenic SNPs，pSNPs）以及基因间SNPs（Intergenic SNPs，iSNPs）等三类。组成DNA的碱基虽然有4种，但SNP一般只有两种碱基组成，所以它是一种二态的标记，即二等位基因（biallelic）。 由于SNP的二态性，非此即彼，在基因组筛选中SNPs往往只需+/-的分析，而不用分析片段的长度。

纯合SNP和杂合SNP是SNP calling软件如GATK或者SAMtools根据测序深度、碱基质量值、比对质量值和基因型质量值等综合判断出来的纯合和杂合，简单来说，纯合SNP可以认为该位点测到的所有reads只是一种碱基类型，杂合SNP为二种或二种以上的碱基类型，不排除特殊位置。

1.6 组装软件

A few aligners are available, including LAST24, BlasR25, BWA-MEM26, GraphMap27,

MECAT28, and minimap229. Only one stand-alone method, PBHoney18, is available to detect all types of SVs from long-read data, although others such as SMRT-SV30 have been proposed for a subset of SV types.

1.7 Sniffles

A particularly innovative feature of Sniffles is its ability to detect nested SVs, such as inverted tandem duplications (INVDUPs) and inversions flanked by indels (INVDELs).

Sniffles using NGMLR with only 10–30× coverage, which recovered around 80% of the calls with precision of ~80% or higher.

1.8 组装过程

基因组组装一般分为三个层次：contig, scaffold和chromosomes。

contig表示从大规模测序得到的短读(reads)中找到的一致性序列。组装的第一步就是从短片段(pair-end)文库中组装出contig。进一步基于不同长度的大片段(mate-pair)文库，将原本孤立的contig按序前后连接，其中会调整contig方向以及contig可能会存在开口(gap,用N表示)，这一步会得到scaffolds。最后基于遗传图谱或光学图谱将scaffold合并调整，形成染色体级别的组装chromosome

image.png

Contig N50：Reads拼接后会获得一些不同长度的Contigs.将所有的Contig长度相加,能获得一个Contig总长度.然后将所有的Contigs按照从长到短进行排序,如获得Contig 1,Contig 2,contig 3…………Contig 25.将Contig按照这个顺序依次相加,当相加的长度达到Contig总长度的一半时,最后一个加上的Contig长度即为Contig N50.

举例：Contig 1+Contig 2+ Contig 3 +Contig 4=Contig总长度*1/2时,Contig 4的长度即为Contig N50.ContigN50可以作为基因组拼接的结果好坏的一个判断标准。

1.9 Scaffold N50

Scaffold N50与Contig N50的定义类似.Contigs拼接组装获得一些不同长度的Scaffolds.将所有的Scaffold长度相加,能获得一个Scaffold总长度.然后将所有的Scaffolds 按照从长到短进行排序,如获得Scaffold 1,Scaffold 2,Scaffold 3…………Scaffold 25.将Scaffold按照这个顺序依次相加,当相加的长度达到Scaffold总长度的一半时,最后一个加上 的Scaffold长度即为Scaffold N50.

举例：Scaffold 1+Scaffold 2+ Scaffold3 +Scaffold 4 +Scaffold 5=Scaffold总长度*1/2时,Scaffold 5的长度即为Scaffold N50.

Scaffold N50可以作为基因组拼接的结果好坏的一个判断标准.

1.10 Biallelic vs Multiallelic sites

A biallelic site is a specific locus in a genome that contains two observed alleles, counting the reference as one, and therefore allowing for one variant allele. In practical terms, this is what you would call a site where, across multiple samples in a cohort, you have evidence for a single non-reference allele. Shown below is a toy example in which the consensus sequence for samples 1-3 have a deletion at position 7. Sample 4 matches the reference. This is considered a biallelic site because there are only two possible alleles-- a deletion, or the reference allele G.

image.png

A multiallelic site is a specific locus in a genome that contains three or more observed alleles, again counting the reference as one, and therefore allowing for two or more variant alleles. This is what you would call a site where, across multiple samples in a cohort, you see evidence for two or more non-reference alleles. Show below is a toy example in which the consensus sequences for samples 1-3 have a deletion or a SNP at the 7th position. Sample 4 matches the reference. This is considered a multiallelic site because there are four possible alleles-- a deletion, the reference allele G, a C (SNP), or a T (SNP). True multiallelic sites are not observed very frequently unless you look at very large cohorts, so they are often taken as a sign of a noisy region where artifacts are likely.

image.png

1.11 Super-enhancer

In genetics, a super-enhancer is a region of the mammalian genome comprising multiple enhancers that is collectively bound by an array of transcription factor proteins to drive transcription of genes involved in cell identity.[1][2][3] Because super-enhancers are frequently identified near genes important for controlling and defining cell identity, they may thus be used to quickly identify key nodes regulating cell identity.

1.12 Enhancers

have several quantifiable traits that have a range of values, and these traits are generally elevated at super-enhancers. Super-enhancers are bound by higher levels of transcription-regulating proteins and are associated with genes that are more highly expressed.[1][5][6][7] Expression of genes associated with super-enhancers is particularly sensitive to perturbations, which may facilitate cell state transitions or explain sensitivity of super-enhancer—associated genes to small molecules that target transcription

every read that spans at least two heterozygous variants can essentially be considered as a 'mini haplotype' that can be assembled into longer haplotype segments by partially overlapping reads spanning the same variable locus。To this end, haplotype-informative reads need to be partitioned into two disjoint sets that represent the two haplotypes。However, all approaches to reconstruct haplotypes from sequencing reads, be it reference-based or reference-free, come with the intrinsic limitation that the distance between subsequent heterozygous markers can be larger than the read length itself. While long-read sequencing (such as PacBio SMRT26 and Oxford NanoPore MinION27), or linked read data (such as those provided by 10X Genomics28) help to mitigate this issue, these technologies still fail to phase over long stretches of homozygosity or centromeres.

1.13 基因的共表达

通常指的是两个基因的表达量呈现相关性。比较常见的例子就是：转录组因子和靶基因间的关系。因为转录因子对它的靶基因有正调控作用，所以转录因子的表达量提高会导致靶基因的表达量也上调，两者往往存在正相关关系。这个正相关关系，可以使用相关系数r^2来度量，这个数值在-1~1之间。总而言之，相关性可以理解为两个元素共同变化，步调一致。

1.14 连锁不平衡（LD）

就是度量两个分子标记的基因型变化是否步调一致，存在相关性的指标。如果两个SNP标记位置相邻，那么在群体中也会呈现基因型步调一致的情况。比如有两个基因座，分别对应A/a和B/b两种等位基因。如果两个基因座是相关的，我们将会看到某些基因型往往共同遗传，即某些单倍型的频率会高于期望值。

例如在群体中（A，a，B，b）各个基因型的频率已知的情况下，各种单倍型的期望频率（AB、Ab、aB、ab）都是可以计算出来。例如，AB的频率=（A的频率）X（B的频率）。但我们实际统计群体中各个单倍型的频率的时候，会观察到某些单倍型的频率会大于期望值，例如下图中的单倍型AB的理论频率是0.12，但观察到的实际频率是0.29。那么说明，基因型A更倾向于基因型B共同遗传。

单体型块 关于染色体位点间相关性的强弱，可以使用LD系数，或者使用一个类似的数值D’来衡量。那些彼此间相关性强的位点，在传代过程中会更少重组，像个“模块”一样倾向于共同遗传。 当然，这种位点间的相关性在染色体上并非是平均分布的。在染色体某些区域，位点间可能会更加紧密相关，共同遗传，这样的区域一般被称为“单体型块”；而某些区域，则位点相关性较弱，重组更加频繁，形成重组区或重组热点。

image.png

1.15 VAF

variant allele frequency 或者 variant allele fraction，对于NGS测序数据来说，就是跟参考基因不同的reads与总的测序reads的比值。正常人的二倍体基因组位点只有杂合或者纯合两种情况，对于纯合那么vaf必然是1，对于杂合，必然是0.5。但是现实测序得到的结果远比这要复杂，尤其是测序深度不够的时候。因为测序本身具有随机性，而且还有很多系统误差。

image.png image.png

在记录SNP的vcf文件中，最后一列GT处只显示了 1/1和0/1这2种类型，说明在call变异位点时，只记录发生变异的位点；对于正常位点0/0不记录；此外，annovar软件注释时也不注释0/0位点。

1.16 StrandPhaseR算法原理

基于两个并行矩阵的二进制排序策略，是一种改进的phase算法，存储从Strand-seq文库单个cell中获得的单倍型信息。输入文件是单个cell的BAM文件；每个Strand-seq文库都定位了单倍型信息的WC区域；不同位置的alleles被识别为W\C reads，以降低单个cell但体型的密度；部分单个cell的但体型被用来生成2个矩阵，矩阵中每一行代表一个cell，每一列代表这个cell中的SNV；2个矩阵分别存放Waston templates和Crick templetes的SNV；矩阵的行按照SNV数量由多到少排列。首先计算每一列的得分（所有变异位点减去最丰富的那个变异），这个分数代表不同cells在指定列的SNV上的差异（类似于每一列进行标准化处理？）；每一列的分数代表了这个矩阵的综合分数，矩阵的综合分数越小，表示所有SNV位点间的一致性越高，即分型越准确。当2个矩阵的分数都计算好后，交换第一行（即第一个cell，含有最多的SNV），这相当于交换了cell的Waston templates和Crick templetes的SNV，以检测是否有更高矩阵得分，即更准确的phased SNV；因此，又重新计算了矩阵得分，如果得分更低则保留之前的得分。这个步骤对所有的cells进行：对两个矩阵中的单细胞单倍型进行排序，以减少每个列中相互冲突的等位基因的数量。

最终的一致单倍型以标准化的VCF格式导出，每个变量位置都有一个指定的Phred质量分数和

熵值，反映了对给定等位基因的可信度。

1.17 What if I find high levels of contamination?

One thing to rule out is sample swaps at the read group level.

Picard’s CrosscheckFingerprints can detect sample-swaps at the read group level and can additionally measure how related two samples are. Because sequencing can involve multiplexing a sample across lanes and regrouping a sample’s multiple read groups, depending on the level of automation in handling these, there is a possibility of including read groups from unrelated samples. The inclusion of such a cross-sample in the tumor sample would be detrimental to a somatic analysis. Without getting into details, the tool allows us to (i) check at the sample level that our tumor and normal are related, as it is imperative they should come from the same individual and (ii) check at the read group level that each of the read group data come from the same individual.

Again, imagine if we mistook the contaminating read group data as some tumor subpopulation! The tutorial normal and tumor samples consist of 16 and 22 read groups respectively, and when we provide these and set EXPECT_ALL_GROUPS_TO_MATCH=true, CrosscheckReadGroupFingerprints (a tool now replaced by CrosscheckFingerprints) informs us All read groups related as expected.

1.18 CNV相关

<u>https://www.jianshu.com/p/eadfc45f1f18?utm_campaign=maleskine&utm_content=note&utm_medium=seo_notes&utm_source=recommendation</u>

TCGA里面主要是通过Affymetrix SNP6.0 array这款芯片来测拷贝数变异！对SNP6.0的拷贝数芯片来说，通常是用PICNIC等软件处理原始数据，就可以得到的segment记录文件，每个样本一个结果，下面是示例结果：

image.png

表明了某条染色体的某个区域内，SNP6.0芯片设计了多少个探针，芯片结果的拷贝数值是多少(这个区域的拷贝数用Segment_Mean)。

通常二倍体的Segment_Mean值为0，可以用-0.2和0.2来作为该区域是否缺失或者扩增，也有人选择0.4作为阈值。

1.19 杂合性缺失

英文是“Loss Of Heterozygosity”，简称“LOH”。

正常情况下，常染色体上的一个区段，都会有来自于父亲、和母亲的各一个拷贝。

当发生杂合性缺失（LOH）的时候，两个染色体上的同一个区段，都是来自于或者父亲、或者母亲的一方，而把另一方的对应区段给丢失了。这就叫杂合性缺失（LOH）。

1.20 单亲二染色体

就是“Uniparental Disomy”，简称“UPD”，是杂合性缺失的一种特殊形式。也就是一对染色体，都是来自于父亲、或者母亲中的一方，而把另一方的对应染色体，全部给缺失了。

1.21 CGH芯片

CGH芯片，也就是“Comparative GenomicHybridization”芯片，“比较基因组杂交”芯片。主要是检测：杂合性缺失（LOH）、单亲二染色体（UPD）、和拷贝数变异（CNV）。

Agilent CGH生物芯片工作的原理：就是把样本的DNA片段化，标上红色荧光素“Cy5”；同时，再把来自几十个正常人的基因组DNA，混成一个标准DNA样本，取同样的DNA量，同样片段化，标上绿色荧素“Cy3”。然后，把这两种标了荧光素的DNA片段，混合在一起，在同一张芯片上进行杂交。接下来进行激光扫描，比较红光荧光与绿光荧光的光强。所得到的光强比值，换算成以2为底的Log值。如果在一个探针上，Log值接近于“0”，也就是说，红光与绿光的荧光光强差不多，那么，可以基本断定，在这个位置，样本中是有2个基因拷贝。如果在一个点上，Log值大约等于“1”，也就是说红光的光强，是绿光的2倍，那么说明，样本在这个位置的拷贝数，可能是标准品的2倍。也就是说，样本在这个位置，可能有4个基因拷贝，比正常情况多出了2个拷贝。同样道理，在一个点上，如果Log值小于等于“-2”，也就是说，红光的强度只有绿光强度的“1/4”，甚至更低，那么说明，样本在这个位置的2个拷贝，可能都丢失了。

1.22 滑窗

对单个样本的bam文件进行200kb的窗口进行滑动计算每个窗口的gc含量（该200kb窗口内每个碱基位点进行统计，再求和），该窗口区域覆盖的reads数量，还有比对的质量值，很容易写脚本进行计算。

GENOME='/home/.../human_g1k_v37.fasta'

bam='ESCC13-T1_recal.bam'

samtools mpileup -f bam |\

perl -alne '{F[1]/200000); F[0]\tGC{F[2]=~/[GC]/;key}+=number{\t}\t}\t}" foreach keys %number}' |\

sort -k1,1 -k 2,2n >T1.windows

得到的结果如下：

image.png

前面两行是窗口的坐标，第几号染色体的第几个窗口，后面3行是数据，分别是每个窗口的测到的碱基数，GC碱基数，测序总深度。

1.23 GC含量与测序深度的相关性

这个是二代测序本身的技术限制，很容易探究到测序深度和GC含量是显著相关的，代码如下：

a=read.table('T1.windows')

a$GC = a[,4]/a[,3]

a$depth = a[,5]/a[,3]

a = a[a$depth<100,]

plot(adepth)

library(ggplot2)

GET EQUATION AND R-SQUARED AS STRING

SOURCE: http://goo.gl/K4yh

lm_eqn <- function(x,y){

m <- lm(y ~ x);

eq <- substitute(italic(y) == a + b %.% italic(x)*","~_italic(r)^2"="~r2,

list(a = format(coef(m)[1], digits = 2),

b = format(coef(m)[2], digits = 2),

r2 = format(summary(m)$r.squared, digits = 3)))

as.character(as.expression(eq));

}

p=ggplot(a,aes(GC,depth)) + geom_point() +

geom_smooth(method='lm',formula=y~x)+

geom_text(x = 0.5, y = 100, label = lm_eqn(adepth), parse = TRUE)

p=p+theme_set(theme_set(theme_bw(base_size=20)))

p=p+theme(text=element_text(face='bold'),

axis.text.x=element_text(angle=30,hjust=1,size =15),

plot.title = element_text(hjust = 0.5) ,

panel.grid = element_blank(),

panel.border = element_blank()

)

print(p)

可以很明显看到GC含量和测序深度是高度相关的：

image.png

1.24 FSHD

DUX4基因位于4号染色体末端的D4Z4区域，该区域由11到100多个重复片段组成，每个重复片段大约3.3 kb长。整个D4Z4区域通常是超甲基化的，这意味着它有大量的methyl groups(由一个碳原子和三个氢原子组成)附着在DNA上。methyl groups的加入会使基因沉默，因此DNA超甲基化区域中激活的基因往往更少。D4Z4区域的每个重复片段都含有DUX4基因，最接近ch4末端的拷贝被称为DUX4，而其他拷贝被描述为"DUX4-like" or DUX4L。

正常人中，D4Z4区域的超甲基化使DUX4-like基因始终保持沉默。DUX4基因在大多数成年人细胞和组织中也是沉默的，尽管它在发育早期和成年男性的睾丸中是活跃的。

DUX4基因(最接近chr4末端的拷贝)位于DNA的调控区域pLAM序列旁边，pLAM序列是产生DUX4蛋白所必需的。chr4的拷贝中，有 functional pLAM的序列被称为4qA序列（此时DUX4基因被激活而表达，使人生病）；而没有 functional pLAM序列的被称为4qB，就不会产生DUX4蛋白。因为每个细胞有2个chr2拷贝，个体可能有2个chr4的4qA序列；2个4qB；或者每个拷贝一个。

（1）面肩肱型肌营养不良症是一种遗传性肌肉疾病，其发病率居于肌肉系统疾病的第三位，发病率为1/20,000。大部分病人在20岁之前已有症状产生，受其影响最严重的是脸、肩、上臂等部位的肌肉，会出现逐渐加重的肌力减退以及肌肉的萎缩。

（2）面肩肱型肌营养不良症分为FSHD1型（占发病人数的95%）和FSHD 2型（占发病人数的5%）。FSHD1型需要由两个遗传因素同时存在才能导致发病：

a. 大约95％的患者（FSHD1型）在染色体4q35区域具有一个称为D4Z4的3.3kb的重复单元的数量缺失，正常人一般有11-150个重复，但是患者只有10个或更少的重复。重复数目与临床表型之间有负相关的关系，重复单元数目越少，发病年龄越早，病情则越严重。

b. D4Z4最后一个重复单元远端需要存在一个称为4qA的变体结构才能导致疾病发生。

1.25 外显子

人类外显子长度平均是200bp（所以默认的bin是267bp，这样可以把比较长的exon给拆分开来）。

1.26 Somatic CNV变异研究方法

随着测序成本的降低以及测序深度的增加，read count 成为最主要的方法。Read count 方法原理是利用一个非重复滑动的窗口去统计覆盖到与该窗口重叠的基因组区域内 reads 数量，从而推断发生 CNV 的位置。 Read count 分析方法包括两个步骤：预处理 (preprocessing) 和分段处理 (segmentation)。 预处理步骤可以对样本比对后的 BAM / SAM / Pileup 等文件进行均一化处理，也可以使用 de-noising 的算法去除 WES 数据中存在的偏好性和背景噪音；分段处理步骤会利用一些统计模型对具有相似 read count 的区域合并去预估CNV的大小，常见的统计模型有circular binary segmentation (CBS), hidden Markovmodel (HMM) 等。

肿瘤样本中 somatic CNV 的检测依然存在一些挑战。基本挑战包括：测序数据质量和测序策略选择。首先，基于read count方法检测CNV，最主要的是寻找基因组某一区段内的reads数量与CNV的关系。然而这种关系会受到样本GC偏好性、数据比对偏好性、实验操作背景噪音以及测序偏好性的影响。其次，在选择WES或TRS测序时，探针的捕获偏好性以及reads在不同外显子区的分布偏好性会影响到CNV检测时的数据统计，会成为CNV检测算法的背景噪音。最重要的挑战是肿瘤样本本身的复杂性，包括肿瘤纯度、倍性以及克隆结构异质性。 CNV 在肿瘤样本基因组上具有广泛性以及多样性，因此 germline CNV 与 somatic CNV 不同之处在于，somatic CNV 可以发生在基因组任何区域且突变频率低。其次，肿瘤样本中污染正常细胞时，会降低 read count 和 read depth 值，使 BAFs 值脱离理论值，影响分段步骤中 CNV 数量估计；基因组非整倍性情况的存在会严重影响肿瘤样本中 BAF 的状态，以及 read count 和 read depth 的基线。最后，肿瘤样本中克隆结构存在异质性，导致一些低频亚克隆结构检出困难，虽然提高测序深度可以帮助检测低频亚克隆，但是想要精准检测 somatic CNV 还是需要综合考虑样本的复杂程度。

Control-FREEC 软件自推出以来人气就很旺，许多高分文章也有引用。该软件可以输入多种类型的预处理文件，并且能够准确区分 somatic 和 germline CNV 变化，尤其适用于肿瘤样本的分析。

1.27 SCNAs

somatic copy-number alterations

1.28 oncogenes and tumor suppressor genes

原癌基因、抑癌基因

1.29 Focal CNVs

(1)Focal CNVs are regions of repeated genetic information that only span a small proportion (<25%) of the chromosome arm (although this does not seem to be a consistent rule), and can contain few genes. CNVs greater than this are termed either 'large-scale' or 'broad'.

(2) focal (length < 98% of a chromosome arm) and arm-level (length > 98% of a chromosome arm)

1.30 Affymetrix SNP 6.0芯片

涵盖超过1,800,000个遗传变异标志物：包括超过906,600个SNP和超过946,000个用于检测拷贝数变化（CNV，Copy Number Variation）的探针，如下：

（1）482,000个SNPs来自于前代产品500K和SNP5.0芯片；

（2）424,000个SNPs的来源包括：HapMap计划中的TagSNPs，X、Y染色体和线粒体上更具代表性的SNPs，来自于重组热点区域的SNPs，以及前代产品未收录的dbSNP数据库中的新SNPs；

（3）202,000个用于检测5,677个已知拷贝数变异区域的探针，这些区域来源于多伦多基因组变异数据库。这些区域可归为3,182个非重叠区域片段，平均每个区域采用61个探针来检测；

（4）744,000个平均分配在整个基因组上的探针，用来发现未知的拷贝数变异区域。

1.31 环状二元分割法(CBS)

全基因组CNVs扫描后面临的一个重要问题就是统计分析．目前已经发展了多种在全基因组水平推算CNVs的软件包和算法模型．其中比较常用的算法是隐马可夫模型(HiddenMarkovModel，HMM)、环状二元分割(Circularbinarysegmentation,CBS)、等级分割(Segmentationalgorithm)、核平滑算法(KernelSmoothingalgorithm)等。

二元分割：在二元分割中，使用 Z-test公式检查出断点所在的位置，然后确定断点所做的分割区域。断点将染色体片段分成2个部分，若断点所分割的2片段的差异度大于片段分割阈值，则肯定有一部分为变异片段。通过循环查找所有显著断点，以至于找到所有CNV变异片段。其基本过程为：在片段中根据Z-test公式查找其最大断点，判断此断点分成的2片段差异度是否大于片段阈值，以确定是否为显著断点。在显著断点分割形成的两侧片段中，再次查找出所有显著断点，最后根据显著断点所划分的变异区域计算其CNV值。

环状二元分割：(Circular Binary Segmentation, CBS)算法是应用在CGH平台上的拷贝数变异检测算法，它是对二元分割方法的修改以提高CNV的检测效果。为检测出隐藏在大片段中的小型变异区域，CBS算法将染色体两端连接起来形成环状分割，通过生成重置排列方法比较非正常数据片段。CBS可以有效地检测出小片异区域。

（1）断点：假设X1,X2,...是一组随机变量，若X1,X2,...XV有符合分布函数 F1,XV,...符合另一个分布函数F2，而F1和F2不同，则索引位置V点即为断点。

（2）BAF：表示基因位点中携带了B等位基因的杂交样本比例。在一个正常样本中，BAF有3种值，即0.0、0.5和1，这分别表示位点的基因型为AA、AB和BB。当计算出BAF的值与这3个值有偏差时，说明该位置可能有拷贝数变异存在。

（3）LRR：是基因位点的实际观测信号强度与期望信号强度比的lbn值。若LRR不等于0，则说明拷贝数有变化。这个函数可以测试出发生变异片段的信号强度和正常强度的偏移量，估计出变异强度。

1.32 肿瘤位点个数

人基因组中平均每1M含有1个突变位点，因此：

WGS中：约3000--4000个肿瘤突变

WES中：约200-300个肿瘤突变

1.33 全****基因组倍增(whole-genome doubling, WGD)

1.34 SINE

散在重复序列是与串联重复序列的组织形式不同的另一类重复序列，是散在方式分布于基因组内的散在重复序列。这类DNA序列一般都是中度重复序列。根据重复序列的长度可以分为短分散重复序列(short interspersed nuclear elements，或short interspersed repeated sequences，SINEs)，，在人基因组中的重复贝数达10万以上。重复序列单元长度在1 000 bp以上的称为长散在重复序列(longinterspersed nuclear elements或long interspered repeated sequences，LINEs)，在人基因组中有上万份拷贝。人类基因组中所有SINE之间的平均距离约为2．2 kb。在结构基因内部，结构基因之间和基因簇内，以及内含子中都有SINEs，但在结构基因的编码区内中还没有发现。

散在重复序列是转座序列，但基因组中多数的SINE和LINE都存在启动子区的缺失，使转座所需的酶无法编码；或者是一端的反向重复序列缺失，使酶无法识别

1.35 共济失调

人体姿势的保持和随意运动的完成，与大脑、基节、小脑、前庭系统、深感觉等有密切的关系。这些系统的损害将导致运动的协调不良、平衡障碍等，这些症状体征称为共济失调。

1.36 junction reads

因为RNA-Seq测序的特性，天然的会有一部分数据延伸到内含子区，这部分跨越外显子和内含子的reads就称为『junction reads』，所以RNA-Seq比对软件需要针对此进行优化，而文章做benchmark也考虑到这点。

1.37 soft clipping

指的是比对的read只有部分匹配到参考序列上，还有部分没有匹配上。也就是一个100bp的read，就匹配上前面20 bp或者是后面20bp，或者是后面20bp比对的效果不太好

2. 染色质免疫共沉淀技术

2.1（Chromatin Immunoprecipitation，ChIP）也称结合位点分析法，研究体内蛋白质与DNA相互作用的一种方法，通常用于转录因子结合位点或组蛋白特异性修饰位点的研究。 将ChIP与第二代测序技术相结合的ChIP-seq技术，能高效的在全基因组范围内检测与组蛋白、转录因子等互作的DNA片段。

image.png

2.2 ChIP-seq技术原理

在生理状态下，把细胞内的DNA与蛋白质交联（Crosslink）后裂解细胞，分离染色体，通过超声或酶处理将染色质随机切割；

利用抗原抗体的特异性识别反应，将与目的蛋白相结合的DNA片段沉淀下来；

再通过反交联（Reverse crosslink）释放结合蛋白的DNA片段；

纯化；

测序获得DNA片段的序列，最后将这些DNA片段比对到对应的参考基因组上。

2.3 CHIP结果之一

[图片上传失败...(image-d23529-1601190421051)]

染色体号

peak起始位点

peak结束位点

peak区域长度

peak的峰值位点（summit position）

peak 峰值的高度（pileup height at peak summit, -log10(pvalue) for the peak summit）

peak的富集倍数（相对于random Poisson distribution with local lambda）

3. 简单理解RNA-seq

3.1 背景

背景1

假设有一群正常的神经细胞（蓝色）和一群变异的神经细胞（红色）

image.png

那么，为什么它们会出现差别呢？是什么遗传机制导致了这个事情呢？

因此，我们需要看一看它们的基因表达差异

我们知道，每个细胞都由一堆染色体组成，每个染色体由一堆基因组成，当然并不是所有的基因都是活跃的，只有一部分基因是可以表达，而表达的中间过程就要经历mRNA转录本，通过高通量测序，我们就能得知：哪些基因是活跃表达的，并且产生了多少转录本（也就是衡量基因表达量的指标）

image.png

背景 3

将正常的细胞测一遍，再将变异的细胞测一遍，得到它们的表达量，我们后来就是比较它们的表达量差异

image.png

可以看出，基因1在两组样本中差异不大或者没有差异；基因2在正常组中基本不表达，而在变异组中表达量很高，二者差别甚大；基因3有差别但比较小

3.2 RNA-seq主要的3步

Step1 构建测序文库

分离RNA =》将RNA打断成小片段 =〉将小RNA片段反转录成DNA =》加接头

接头两个作用：测序仪识别；允许一台测序仪同时运行多个样本，提高性价比

但是需要注意：加接头的过程是随机的，并不是所有的接头都被加上，有些反转录的DNA片段没有加上接头

=》PCR扩增（只有加上接头的测序片段才能被扩增）=〉质量检查QC（看下文库的浓度和片段长度）

image.png

对文库进行测序

一块测序板上（想象下载玻片，其实人家真名是Flowcell）能包含多于400,000,000个片段，垂直于测序板排列。

测序仪有四种颜色的荧光探针A、T、C、G，与测序片段上碱基互补，结合上就“放烟花”表示庆祝

image.png

测序就是这样，结果就得到了raw data，就是fastq数据

Step 2 原始数据处理

质控=》过滤garbage reads=〉比对到参考基因组=》再数一下每个基因比对上多少reads

garbage reads：

有些时候接头并没有加到测序片段，而是他们直接结合，也能进行测序，但测得结果是没用的

比对到参考基因组

先将大的基因组序列打断成许多小片段，然后为了方便接下来寻找这些片段，需要对他们进行构建索引index（目的就是标注每个小片段的位置）

再将测序的reads和基因组一样，也是打断成小片段，然后把它的小片段比对到基因组的小片段上，比对上的会给出位置信息

image.png

统计reads数得到表达矩阵

就想这样：第一列是基因名（人类基因组有大概2w基因，因此大概有2w行）

其他列是每个测序样本比对上的数量（6-成百上千不等），这里的6的考虑的是处理对照各3个重复，即Bulk-seq；大样本量的RNA-seq比如Single-cell，每个细胞都是一个样本，因此成百上千

每一行都是原始的统计值，每个基因在每个样本中被抓到多少次

image.png

标准化表达矩阵

进行标准化的原因是：某些样本可能本身测序质量就差，但并不代表人家没东西；或者有的样本测序的时候加的浓度比较高，因此统计时占优势，但并不公平！

因此需要让大家在同一起跑线

image.png

Step 3 可视化

比如PCA分析，看看样本之间能否区分开，另外可以排除明显不对的样本，比如这里的wt2

image.png

然后看差异表达基因（就是正常与变异样本的差异）

红色是差异的，黑色是共同的

image.png

如果发现了感兴趣的差异基因，怎么办？

这个基因是你研究的，接下来通过实验验证

对这个基因不熟悉，只是感兴趣，就可以做GO、KEGG注释，看看它在正常还是变异样本中有富集。

4. 二代测序数据拼接之原理篇

4.1 拼接基本原理

拼接可以分为基因组和转录组拼接，基因组拼接对数据量和测序深度要求更高，而转录组用平时的 RNA-seq 数据就可以。目前手里更多的是 RNA-seq 数据，所以做的也是和转录组拼接相关的内容。

无论拼基因组还是转录组，归根结底都是拼 DNA 序列。拼接最简单的逻辑就是一个由让由长变短再由短到长的过程，基因打断成一个个片段进行测序，生成测序数据，然后再用reads拼成contigs再拼成scaffolds。

从 reads 到 contigs 的过程中，需要进行多序列比对并将一致的 reads（consensus sequence）拼接起来来生成contigs，其中最大的问题是基因组上存在大量重复区域，会对拼接带来困扰。Scaffolds数据则是通过pair end reads 信息来判断contigs 的顺序、方向和相邻 contigs 之间的缺口（gap）大小来生成的。从contigs到scaffolds是一个排序和定向的过程。

4.2 拼接常用算法

目前常用的拼接算法都是基于数学中的图论思想（Graph theory）产生,其中图论中的两个点表示两个read，而两点之间的连线表示两条read的重叠区域。拼接要做的事情就是在所有的路线中找最优解，类似于小时候玩过的一笔画问题。

image.png image.png image.png

如上图左所示，一个简易的基因组产生了ABCD……若干read，在理想的情况下我们可以根据所有reads彼此的重叠区域重构出圆圈所表示的基因组。如果简化成图来表示则应该是上图右的所示，所有黑连接的是正确的基因组，但实际情况是基因组有很多区域比较相似，以至于序列间会产生本没有的联系（如红线所示）。

在进行上述分析过程时，需要把所有的reads都进行比对，以便找到重叠区域，这个步骤非常耗费计算资源。Graph画出来之后的问题就是如何从中得到最优路径，即从有多种reads组合方式中找到从合适的一个形成contig。

下面是寻找最优路径的常用算法。

Greedy extension

贪婪算法（贪心图）是早期提出的拼接算法。首先选定初始read, 然后找和其重叠区域最高的read进行延伸，直到拼接后的read两端都不能再进行扩展为止。

每一次都是从最优匹配开始，然后次优匹配，到不能匹配时停止。这样一来，贪婪算法通常会得到局部最优解，而不是全局最优解。因此，这种算法在遇到重复序列时会出现比较大的问题。

A greedy assembler compares all pairs of read fragments continually and replaces any pair with sufficient overlap between the edges with a combined sequence. The assembly completes, when the edges of the remaining sequences do not have any significant overlap.

Overlap Layout Consensus

OLC图算法主要是用来针对长reads序列拼接，如一代测序数据（三代测序数据），简单理解就是把测序产生的长序列用彼此之间的overlap区域连接起来。对于数据量很大的数据或者全基因组数据来说，形成的olc图非常复杂，会消耗大量内存。

OLC算法共有三步：

Overlap

对所有reads计算任意两条之间的重叠区域，挑选出满足筛选条件的reads。这里区别与贪婪算法，会先把所有overlap都找到。这一步，通常会将一个reads分成若干个长度比较短的序列（kmer/seed/word），要求是每个片段序列之间至少有若干个碱基的重叠区域。

layout garph

简单化过程。对reads进行排序，确定reads之间的位置，建立overlap图，将重叠的reads组合成contig。

Consensus

在已经建好overlap图的基础上，将所有的read序列排列起来，找一条从起始节点到终止节点的最佳近似路径使得最终路径将会遍历一次重叠区域中的每个节点，相当于对初始的reads集合中全部序列进行重构得到目标基因序列。

de Bruijn graph

De Bruijn 图是目前最常用的二代测序拼接算法。比较流行的拼接软件如 Velvet、Abyss 和 SOAP denovo 都使用该算法。

与OLC不同之处在于，这个算法将已经非常短的reads再分割成更多个kmer短序列（k 小于reads 序列的长度），相邻的kmers序列通过（k-1）个碱基连接到一起（即每次只移动一个位置），进而降低算法计算重叠区域的复杂度，降低内存消耗。

这里的kmer首先不能太短，比如2个碱基肯定拼不出来基因组，它的长度既需要能够使其携带足够的基因组的信息，也要短到可以进行后续的错误矫正。除此之外，一个read中小的片段被分割之后还不会丢失原来reads 的前后位置信息。

总体而言，该算法将reads打断成长度为K的核酸片段，再用Kmer间的overlap关系构建DBG，最后通过DBG得到基因组序列。

拼接步骤通常包括：

§

构建DBG图，将read分割为一系列连续kmer

§

如下图，1和2两个序列会拼接成一条长的序列，在引入第三条序列后会出现一个圆环。

image.png

下图是一个最简单的DBG拓扑结构，两球一线（一进一出）代表相邻的两个kmer，圆圈则代表有多种连接方式。

image.png

下图是三种常见的DBG结构，分别代表了拼接过程中的不同情况。

image.png

§

合并DBG图

§

合并路径中出度入度唯一（one incoming and one outcoming ）的节点，去除段末端，低覆盖度节点和泡状结构。

image.png

§

构建contig

§

寻找最优路径（经过每个节点且仅经过一次），最优路径对应的碱基序列构成一个contig

§

构建scaffold

§

通过PE reads 位置信息确定contig之间的相对位置和方向，组装contig，填充contig之间的gap，得到scaffold序列。

两个注意事项

（1）当把双链信息考虑进来之后，可能的连接情况就会增加，因为任何一个节点都可能和某个节点的反向互补序列相连。

（2）有重复区域时的情况如下图所示，可以发现，一个确定的genome只可能有一个DBG，但反过来一个DBG不一定找到的是唯一的潜在基因组。

image.png

kmer和内存

在拼接相关的文章中，kmer是出现频率非常高的一个词。而kmer在整个生物信息分析过程中的用处也是非常之多。

k值越大可辨别更多的小重复序列，越容易把DBG转换为唯一的序列，但得到的拼接过程含有更多的gaps；小的k值对应的DBG能够得到较好的连通性，但是算法的复杂度会提高，repeats序列处理会更复杂，增加了错拼的可能性。

在拼接数据预处理软件khmer的文献中有这样一段关于kmer和内存大小与处理结果关系的描述：

The interaction between these three parameters and the filtering process is complex and depends on the data set being processed, but higher coverage levels and longer k-mer sizes result in less data being removed. Lower memory allocation increases the rate at which reads are removed due to erroneous estimates of their abundance, but this process is very robust in practice

kmer越大需要的内存就越多，所以计算机的内存大小也会限制kmer的取值。这里需要说明的是，输入数据的多少不会影响memory用量，但是输入数据的错误越少，占用的内存也就越少，假设所有测序数据都没有任何错误，那么DBG的大小并不会因为测序深度的增加，因为不需要将因为几个碱基不一致的kmer存入到DBG中（下一部分会具体提到）。至于需要多大的RAM则取决于DBG的大小和组装基因组的大小。

另外，在拼接的过程中尽量避免使用偶数kmer，否则容易是kmer产生回文序列，特别是在链特意性的数据中。

在平时分析中，一般会设置一个kmer的梯度（21，23，25，27,2931），来解决DBG算法loss of read coherence的问题。然后从中选择最好的结果。另外，还有一种说法是在进行拼接过程时，kmer应该选择read长度的1/2到2/3大小，否则可能拼接出过多的Contig。这一点，也可能是我们平常使用trinity拼接时拼出Contig 过多的原因，trinity的默认拼接大小是25。上限是32？（有待确定）。如果kmer有上线，是否也可以考虑在预处理的时候，处理的力度大一点，把序列截短一些？

拼接的干扰因素

在实际情况中，拼接往往是在覆盖度不均匀且含噪声的数据中进行，这为拼接带来了三个方面的困难：

增加了大量假kmer从而提高了对缓存的要求

错误的reads通过增加tips，bubbles和corss-links等改变了DBG的结构

不统一的read覆盖度使得拼接参数对拼接结果有非常大的影响

假kmer

在一次测序得到的数据中，kmer matches 的数量和测序深度以及read长度相关。假设在完全没有测序错误的数据中，read长度是100，测序深度是50X，选取kmer值为21，那么一个只匹配基因组一次的kmer 出现的次数应该是（基因组的每个位置被测了50次，100bp的read有80个21bpkmer）,如果匹配基因组两次的kmer应该出现80次。因此，峰值应该在40，80的位置有一个小峰。

但是，当存在测序错误的时候，会在匹配1次的位置出现大量的kmer,就是由于测序的误差导致的。为什么说是测序错误呢，因为在50x的测序中只出现了一次，如果一个read中有一个碱基错了，那么这个read就会产生21个错误的21kmer。更大的问题是，随着测序深度的增加，这样错的kmer数量也会增加，

改变结构

错误的reads 会为DBG引入三种类型的错误

tips

所谓tips指一个小分支，下图所示，我们有5条10bp的reads，其中第5条有一个碱基测序错误。使用7kmer会产生11个节点。如果用前4个read来拼是正常的，一旦引入第五条就会因为一个错误的碱基出现一条错误的有三个节点的分支。

image.png image.png

在大量拼接的过程中，大量的read会匹配到正确的位置，一小部分会比配到错误的位置，因此可以对错误的tip进行清除。

image.png

bubbles

当read和kmer相比足够长时，错误可能出现read中间。此时会出现bubble的情况

image.png

需要注意的是，除了测序错误以外，可变剪切和snp以及插入缺失等也会导致tips的出现，因此增加了拼接的难度。因此，在进行拼接之间有必要对原始数据进行预处理。

覆盖度不均一

在实际拼接过程中，会去除一些低频率的kmer，这一操作在删除了大量错误kmer的同时，我们也不可避免的删除了很多是因为覆盖度低所以出现次数低的kmer。

综上，如果调低cutoff，高覆盖区域出错的可能就高，但是低覆盖区域的质量提升。kmer长度增加会使低覆盖区域更加分散，因为kmer的覆盖度会因为低于设置的cutoff被删除。

image.png

随着kmer的增加，分支会逐渐减少，DBG会越来越趋向于线性。对于一个很大的kmer,可能就完全线性，同时按照染色体分开。拼接质量通常会随着k值的增加先变好然后再变坏，因为这个过程中存在两种竞争性的过程。一方面，kmer的增加可以更好的处理重复区域，但是另一方面，由于覆盖度的原因，kmer的增加在一些区域使得他们出现的次数越来越低直到低于筛选的阈值。所以，kmer 的选择非常重要。

image.png

基因组拼接和转录组拼接

image.png

不同的拼接内容需要不同的拼接策略，其原因如上图所示，即不同的数据产生的DBG结构和覆盖度不同。

对于转录组拼接来说，如果假设一个转录组只有两个基因而且这两个基因没有重复区域，那只需要构建两个没有关联的DBG就可以了。

image.png

对于实际数据来说，一个转录组包含了成千上万个基因。且大多数基因没有overlap，我们构建的DBG实际上是众多不相关grphs的集合，每一个图都代表一个基因或者一个基因家族。下面的示意图显示构建了五个基因的图并且也展示了其覆盖度的差别，其中基因2和3可能是两个可变剪切或者一个家族非常类似的两个基因。

image.png

通过上图也可以发现，转录组拼接最大的问题在于kmer数目的高低很大程度上取决于某个基因的表达量高低。如果把kmer cutoff 设置的高一些，那么低表达的基因就可能拼不起来。

在基因组拼接过程中，kmer coverage 通常是单峰或者是双峰，但是转录组拼接不同。转录组拼接中，每一个基因的kmer分布可能是一个单峰，peak的位置取决于基因表达量的多少。而整体的分布则是所有这些单基因分布的总和。

转录组拼接时应该注意这样量个问题。首先，转录组拼接已经不满足DBG算法中覆盖度均一的假设，表达量非常高的基因可能主宰拼接的结果，因此需要能调整这种覆盖度的差异。另外，转录组拼接不用太担心低复杂度区域，不会有很多重复区域出现。

转录组拼接常用软件：trinity；SOAPdenovo-trans(华大)；Trans-ABySS；bridger

拼接质量的评估

基因组层面的拼接质量，一般会比较看重长度相关的指标。比如最大长度、平均长度、拼接后的总长度和 contig N50长度。

Contig N50 指 reads 拼接后获得一些不同长度的 contigs，将所有的 contigs 长度相加得到总长度。然后将所有的 contigs 按照长度从长到短进行排序，再将 contigs 按照这个顺序相加，当长度等于 contigs 总长度的一半时，最后一个加上的 contig 的长度称为 contig N50。对于总长度不同的两个拼接数据，直接对比N50 的数值没有什么意义。

对于转录组拼接而言，并不是越长越好，我们更在意的是拼接的质量，方向和回帖率等等信息。如果我们在转录组拼接过程中使用了kmer=25这个参数，在拼接好后应该用拼接的fastq文件mapping回拼接好的转录组，测试mapping效率，这里推荐使用salmon软件，需要注意salmon中的kmer应该和拼接时采用的kmer保持一致。

另外，transrate是一个专门用来评价拼接质量的软件。在后续的实际应用部分会有介绍。