结构变异（Structural Variation）检测方法探析

结构变异一般是指染色体重组，在基因组上距离很远的两个基因发生融合，形成了新的编码序列。

这样的结构变异往往会导致比较严重的疾病，所以通过NGS的方法检测结构变异对于医学研究有很大用途。

最常见的思路是，对于已知的基因之间的融合，比如geneA和geneB， 分别进行比对，然后对生成的sam文件进行解析。

这里选用NextGenMap进行比对，方法如下：

ngm -r gene_A.fa -1 ${fq1} -2 ${fq2} -o geneA_region.sam -t 8
ngm -r gene_B.fa -1 ${fq1} -2 ${fq2} -o geneB_region.sam -t 8

-- 采用Python形式参数优化SAM解析过程

python中形式参数可以使用户解析多个变量，示例如下：

def p(*args):
    for i in args:
        print(i)

p(1,2,3)

这里要解析的是对geneA和geneB 比对产生的SAM，然后找出同时比对到两个基因上的reads， 根据reads split的位置判断是否是融合基因。

#!/usr/bin/env python
# -*- coding: UTF-8 -*-

from collections import defaultdict
geneA_sam= './geneA_region.sam'
geneB_sam= './geneB_region.sam'

def getReadsN(*files):
    d = defaultdict(dict)
    for file1 in files:
        opFile = open(file1, 'r')
        for line in opFile:
            if not line.startswith('@'):
                xl = line.rstrip().split('\t')
                if xl[5] != '*':
                    chrName = xl[2]
                    d[xl[0]][chrName] = d[xl[0]].get(chrName, 0) + 1
    
        opFile.close()
    return d
    
testD  = getReadsN(geneA_sam, geneB_sam)

geneNames= ['X1' , 
            'C1' ,
            ]
outFile = open('Reads_ID_stats.txt', 'w')
outFile.write('Reads-ID\t' + '\t'.join(geneNames) + '\n')
for i in testD:
    if len(testD[i]) >  1:
        outFile.write(i)
        for j in geneNames:
            
            outFile.write( '\t{0}'.format( testD[i].get(j, 0) ) )
        outFile.write('\n')
outFile.close()

在结果中 会列出比对到两个基因上的Reads ID，后续的分析需要去找出比对的CIGAR上的位置

比如，Reads1 在geneA的比对信息是39M111S， 在geneB的比对信息是110M40S， 那么Reads1 很可能就代表geneA和geneB的一个融合位置。

M的含义是Match， S的含义是SoftClipping

CIGAR的含义在SAM文件的定义中有解释，参考 https://samtools.github.io/hts-specs/SAMv1.pdf

如下：

Op   BAM    Description    Consumes query    Consumes reference

M   0   alignment match (can be a sequence match or mismatch) yes yes
I      1   insertion to the reference yes no
D    2   deletion from the reference no yes
N    3   skipped region from the reference no yes
S    4   soft clipping (clipped sequences present in SEQ) yes no
H    5   hard clipping (clipped sequences NOT present in SEQ) no no
P    6   padding (silent deletion from padded reference) no no
=    7   sequence match yes yes
X    8   sequence mismatch yes yes

比较常见的就是M, S, I, D，其中I代表Insertion， D代表Deletion

目前检测SV已有很多发表的软件，比如novoBreak， 在 ICGC-TCGA DREAM 8.5 Somatic Mutation Calling Challenge 上得过奖。

还有GROC-SVs, 适用于10X genomics产生的长片段Reads。

以下是结构变异的示意图

Figure 1 表示的是常见的Break End 类型的结构变异

Figure 2 表示的是Chr2 的 321682 position与Chr13 的 123456 position 之间产生了有插入序列的结构变异，这种结果用本文提供的结构变异检测思路就难以检出，需要设计更好的算法。

在记录变异的VCF文件格式中，已有各种类型的SV的记录方法， 参考  http://samtools.github.io/hts-specs/VCFv4.2.pdf

 https://support.10xgenomics.com/genome-exome/software/pipelines/latest/output/sv-vcf