外显子数据分析之基因融合factora
作者，Evil Genius~~~
美国人不按照套路出牌啊，外显子CNV分析原理及各种注意事项刚给下面的人培训完，结果今天考核的基因融合，真是出事了。
让我想起了刚入职参加工作的时候，也是培训天天考核，考核是随机并且没有目的性的，所有内容都在考核范围，有一次我也考核没过，心里那种害怕还记忆犹新。只不过现在站在台上的是我，培训台下一群5年前的“我”。
关于外显子基因融合一般采用factora或者genefuse，我们这里采用factora，感兴趣的可以查阅genefuse。factora文章在FACTERA: a practical method for the discovery of genomic rearrangements at breakpoint resolution
基因融合常见的三种发生机制：
1）Chromosomal Translocation，染色体易位。
2）Interstitial deletion，中间缺失。
3）Chromosomal Inversion，染色体倒位。
知识背景
- 一般来说，基因融合是指基因组层面的融合。但转录组层面也可能发生融合，主要是由于两个不同基因转录产生的RNA，由于某种原因融合在了一起，形成新的融合RNA，该RNA可能编码蛋白，也可能为非编码。而基因组层面产生的融合基因，根据融合的情况，可能表达，也可能不表达（如破坏了启动子区域或其他原因）。
- 导致基因融合的DNA重排是人类恶性肿瘤中获得性结构变异的主要原因。
- 许多结构重排是致癌驱动突变，并且越来越具有治疗靶向性
- 由于其独特的连接序列，融合产物也可以作为cell-free DNA中肿瘤负荷的非常敏感的生物标志物，这些cell-free DNA不断地脱落到各种体液中。
- FACTERA，一个软件工具，用于发现基因组重排，包括易位，倒位和缺失。
METHODS
- 作为输入，FACTERA需要(i)由能够“soft clipping”的比对工具绘制的paired reads的BAM文件（BWA）; (ii)基因组坐标，用于通过基因、外显子或其他基因组信息的位置控制融合发现的分辨率;(iii) hg19.2bit参考基因组，以实现快速序列检索。
①序列比对软件（如BWA）在将序列比对至参考基因组时，会对异常序列进行处理，如果某段序列左右两端能同时匹配到基因组的不同位置，软件会对序列进行切割。通过遍历bam文件中的异常序列，识别出一些有着明显断层的异常序列簇。
②找到这些异常序列簇并对异常序列的左右两段序列分别进行处理，左端序列以10碱基为滑窗，并对滑窗内序列做K-mer（K = 4）处理，结果存入HASH表中，右端序列以4碱基为滑窗，得到短序列集合，若HASH表和短序列集合之间是否存在相似性，那么这就是一个潜在的融合位点。
③序列中出现的重复碱基为断点的判定增加了一定的难度，FACTERA默认选择最左侧的位点作为断点。考虑到反向互补序列的存在，融合基因序列在识别过程中，会出现4种方向，但最终呈现时，只有两种方向（GENE1-GENE2，GENE2-GENE1）。在对断点校正之后，FACTERA会使用BLASTN将异常序列簇以及比对至参考基因组失败的序列与生成的融合基因序列进行比对，从而验证该融合位点的准确性。 - 截断(或“soft clipping”)reads的clipped boundaries代表潜在的融合断点。
- R1和R2有四种可能的取向，然而，只有上图D所示的1a和2a病例可以产生有效的融合，因为它们的reads具有面向相反方向的soft-clipped sequences。
- 最后，为了在read comparison和断点调整后验证候选融合，FACTERA使用BLASTN将所有 soft-clipped a和未映射的reads与每个候选融合序列(在断点上下填充500 bp)对齐。保留至少具有95%标识且超过90% input reads length(默认情况下)的reads，and reads that span or flank the breakpoint are enumerated。
-
除了基本算法之外，还实现了几种算法来提高性能。
融合的表示方法
一、断点方式（SVTYPE=BND）
在描述SV的VCF规范中，其中之一就是详细的描述SV的断点，并在INFO列标明SVTYPE=BND，这时，SV的具体发生情况就有四种类型：
如上图所示，VCF的BND表示法遵守以下的规则：
- 首先列出染色体和第一个断点的位置，第二个断点的染色体和位置在ALT列标出（如上图的1:800）
- 第二个断点的方向由
[的方向决定，[1:800[表示从给定位置向右延伸，]1:800]表示从左边开始延伸并在指定位置结束 - REF 列
N相对于[1:800[和]1:800]的位置即为第一个断点1:500相对于第二个断点的位置 -
N可能是某一个特定序列，这取决于REF列
符号表示法的易位（<TRA>）
VCF中描述SV的另一种方式是“符号表示法”。在VCF的ALT列中相应的标签，如<INV>表示倒位，<DUP>表示重复。
由于易位相对复杂，因此这里的<TRA>表示的也比较模糊，仅有标签是不够的。所以就有了在INFO列增加相应的标签（CHR2表示第二个断点的染色体，END表示具体位置）表示第二个断点的位置和方向。最后还需要的一个信息就是两个断点间的相邻关系，这里用CT表示。
所以，两种表示方法间的对应关系就有了：
| BND | <TRA> with CT INFO field |
|---|---|
| 1 500 . N N[1:800[ | 1 500 . N <TRA> ... CHR2=1;END=800;CT='3to5' |
| 1 500 . N ]1:800]N | 1 500 . N <TRA> ... CHR2=1;END=800;CT='5to3' |
| 1 500 . N [1:800[N | 1 500 . N <TRA> ... CHR2=1;END=800;CT='5to5' |
| 1 500 . N N]1:800] | 1 500 . N <TRA> ... CHR2=1;END=800;CT='3to3' |
其他类型的符号标签（<DEL>, <INV>, <DUP>）
下图展示了删除，倒位和重复的几种情况：
由此可见，删除和重复只有一种相邻关系，而易位则会存在两种情况，两种表示法的对应关系如下：
| Symbolic Call | As BND call(s) |
|---|---|
| 1 10 . N <DEL> ... END=20; | 1 10 . N N[1:21[ |
| 1 10 . N <INV> ... END=20; | 1 10 . N N]1:20] |
| 1 1 . N <DUP> ... END=10; | 1 1 . N ]1:10]N |
在合并多个样本或工具生成的SV VCF时，mergevcf 和 SURVIVOR 均会将符号表示法转为BND表示法，因此，很有必要搞清楚它们之间的对应关系。
基因融合鉴定的一些常见术语(A) Intact exon (IE) type andbroken exon (BE) type fusion transcripts; (B) spanning read, split readand anchor length; (C) short and long insert size of DNA fragment forsequencing.
这些常见的术语有：
-
1）Intact exon （IE） type fusion，是指融合后完整的保留了原来的外显子，未影响原来的外显子结构。如上图A中Gene A的Exon2和Gene B的Exon1融合后完整的保留了两个外显子的序列；
-
2）Broken exon (BE) type fusion，是指融合后没有保留原来完整的外显子序列。如上图A中 Gene A的Exon3的部分序列和Gene B的Exon2融合在一起，融合后的新基因中，来自Gene A的Exon3丢失了部分序列；
-
3）Breakpoint，是指两个融合基因在基因组上发生融合的位置，如上图B中Gene A（蓝色）和Gene B（绿色）融合的位点；
-
4）Spanning read，是指跨越融合位点分别匹配到两个融合基因的paired-end read，比如上图B中的匹配到Gene A（蓝色）和Gene B（绿色）的一对read；
-
5）Split read，是指恰好匹配到融合位点上的read，具体如上图B中右侧图所示；
-
6）Anchor length，是指跨越融合位点的read左端和右端的长度，具体如上图B中右侧图所示；
-
7）short insert size，一般是指双端测序paired-end sequencing中，两个read中间间隔的较短距离，一般为几百bp；
-
8）long insert size，一般是指双端测序mate-pair sequencing中，两个read中间间隔的较长距离，一般为几kb甚至更长；
FACTERA的输出文件(这块儿不翻译了，翻译了感觉表示不出来真实的意思)
Each FACTERA run produces 9 main output files, each of which is described below:
-
parameters.txt = all input files and parameter values.
-
discordantpair.depth.txt = ranked list of discordant read clusters.
-
disordantpair.details.txt = discordant read positions.
-
fusiontargets.bed = bed coordinates of candidate fusions – used to restrict search space for soft-clipped reads.
-
blastreads.fa = used to build blast database of soft-clipped, improperly paired, and unmapped reads.
-
blastquery.fa = file used to search individual candidate fusion sequences (query) for hits in blastreads.fa (target database).
-
fusionseqs.fa = all detected breakpoints with 500bp of additional flanking sequences.
-
fusions.bed = bed output for detected fusions. Useful for comparing runs or somatic vs germline (column 4 is fusion ID).
-
fusions.txt = all detected fusion events, including details, described below:
| Field | Description |
|---|---|
| Est_Type | Estimated structural variant type: TRA = translocation; INV = inversion; DEL = deletion; '-' = not determined |
| Region1 | Name of genomic region closest to breakpoint 1 (e.g., gene 1, exon 1, etc.) |
| Region2 | Name of genomic region closest to breakpoint 2 (e.g., gene 2, exon 2, etc.) |
| Break1 | Chromosomal breakpoint 1 |
| Break2 | Chromosomal breakpoint 2 |
| Break_support1 | Number of reads supporting breakpoint 1 |
| Break_support2 | Number of reads supporting breakpoint 2 |
| Break_Offset | Breakpoint adjustment in bases (e.g., owing to microhomology) |
| Order1 | Orientation of read clipping with respect to breakpoint 1: CN, clipped followed by not clipped; NC, vice versa |
| Order2 | Same as Order1, but for breakpoint 2 |
| Break_depth | Number of breakpoint-spanning reads |
| Proper_pair_support | Number of properly paired and previously soft-clipped reads that map to fusion |
| Unmapped_support | Number of previously unmapped reads that map to fusion |
| Improper_pair_support | Number of previously discordantly paired reads that map to fusion |
| Paired_end_depth | Total number of paired-end reads that flank breakpoint |
| Total_depth | Mean total depth for regions flanking both breakpoints (+/-500bp by default) |
| Fusion_seq | Estimated fusion sequence (50 bases flanking breakpoint by default) |
| Non-templated_seq | Non-templated (i.e., non-reference) sequence segment (if any) enclosed in brackets |
