在做转基因动植物时，相信许多实验室的小伙伴们都遇到过＂我转成功了吗？＂我转到哪儿去了？＂这样的困惑，那么，你是否还在传统的酶切，PCR扩增，构建载体，然后一代测序等分子实验技术来做验证？

你是否知道，利用重测序的方式检测外源基因的插入与否以及确定外源基因插入位点则是一个非常不错的选择呢。来来来，听小编给你娓娓道来。。。

01

外源基因插入位点检测原理

将测序reads比对到参考基因组和外源序列，根据比对结果文件找出下列两类 reads：

第一类：一端reads比对上参考基因组序列，另一端reads比对上外源序列；

第二类：两端中任何一端reads一部分序列比对上参考基因组序列，另一部分比对上外源序列。

检测原理示意图如下：

图1 橘色部分为外源序列，绿色部分为参考基因组。需要从比对结果中找出上面的两类reads，其中第二类reads可以准确定位插入位点。

02

外源序列插入位点的检测步骤

１）对下机原始数据进行过滤得到clean reads；

２）将外源序列与参考基因组序列进行合并得到ref.genome.fa文件；

３）Clean reads与上述fa文件进行比对得到sam文件；

４）将外源序列与参考基因组进行比对，排除同源性的影响；

５）评估外源序列的测序深度和覆盖度；

６）从bam文件中筛选出比对到外源序列的reads，进行组装；

７）组装序列与外源序列进行比对，评估组装结果；

８）组装序列与参考基因组进行比对，确认插入位点。

03

具体步骤介绍

1）下机原始数据会包含adapter，一些低质量的reads以及adapter污染的reads，一般使用fastqc进行质控，利用cutadapt（http://cutadapt.readthedocs.io/en/stable/）进行过滤，也可以用trimmomatic（http://www.usadellab.org/cms/?page=trimmomatic）进行过滤，如有小伙伴拿到的是原始数据可以使用这些软件进行过滤，一般使用默认参数即可,trimmomatic命令

java -jar trimmomatic-0.35.jar PE -phred33 input_forward.fq.gz input_reverse.fq.gz output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

2）在测序界来说，几乎所有后续分析都要基于序列比对，寻找外源序列插入位点也是如此，利用bwa（http://bio-bwa.sourceforge.net/bwa.shtml）和samtools（http://www.htslib.org/doc/samtools.html）获得bam文件

bwa index –a is ref.genome.fa#is 是bwa默认的算法，当database大于2GB时不可用
bwa index –a bwtsw ref.genome.fa#当database大于2GB可选择此算法
bwa mem –t 4 ref.genome.fa example.1.fq example2.fq > example.sam # 生成sam文件，-t 表示线程数
samtools view -bS example.sam > example.bam #sam格式转bam格式

3）对于第4步主要是评估外源序列与基因组之间的同源性（假如同源性很高，怎么判断测序reads是来自外源序列还是参考基因组？），blast（https://www.ncbi.nlm.nih.gov/books/NBK52640/）软件即可

formatdb -i ref.genome.fa -p F #利用参考基因组构建数据库,-p:建库类型
blastall -p blstn -i -d ref.genome.fa -m 8 -o example.out #-m表示输出格式

一般选择m8表格形式，还有其它格式如m7的xml，简单看一下m8输出格式：

图2 各列意义：query名/subject名/identify/比对长度/错配数/空位数/query比对起始坐标/query比对终止坐标/subject比对起始坐标/subject比对终止坐标/期望值/比对得分

4）第五步主要是检测有无reads比对到外源序列，假如没有reads比对到外源序列，材料就不是真正的转基因材料； 如有reads比对到外源序列,需进一步统计外源序列覆盖度与覆盖深度，评估结果可靠性，参考如下：

samtools mpileup -f ref.genome.fa example.bam >example.txt # -f参考文件

简单看一下example.txt格式：

图3 各列意义：参考序列名/位置/参考碱基/比对上的reads数/比对碱基/碱基质量/

根据example.txt第2列和第4列就可以很容易得到外源序列的覆盖度与覆盖深度了。

5）讲第六步之前，大家需要了解bam文件格式，（http://samtools.github.io/hts-specs/SAMv1.pdf），了解格式之后就可以很轻松的把比对到外源序列的reads 信息提取出来，之后就可以去clean reads中提取比对到外源序列的reads了。

samtools view example.bam scaffold_1 > scaffold1.sam #此处假设外源序列名字是scaffold_1

6）接下来需要对筛选的reads序列进行拼接组装，可选择的软件也有许多，velvet，soap denovo ，spades（http://cab.spbu.ru/files/release3.11.1/manual.html#sec3.2

）等，简单点说就是将短reads拼接成更长的contigs，填补gap，由contigs再到scaffolds，reads拼接完成以后,由检测步骤中的7,8（同样使用blast软件）,就可以轻松得到外源基因的插入位点了。

Spades.py -1 example1.fq -2 example2.fq -o d/ # -1 上游reads -2 下游reads -o 输出目录

7）最后一般会使用igv进行截图验证,先随便来一张,如果在下图中能看到箭头所示reads,那就要注意了,很可能那就是你要的结果。

图4 igv截图验证

如果我们对基因组某些区域感兴趣,需要查看这些区域的reads覆盖情况,当然不可能一张张去手动截图,只需要编辑一个igv.batch格式的脚本即可

snapshotDirectory #截图需要保存的路径，
goto Chr01:41,804,384-41,816,127 #需要查看的基因组位置
snapshot 1.png #生成图片的名称
goto Chr01:41,832,182-41,843,925
snapshot 2.png
goto Chr01:41,877,711-41,889,454
Snapshot

然后IGV工具栏Tools --> Run Batch Script,就可以轻松批量截取。如果基因组比较大，查看区域也比较多，不太可能手动编辑如下脚本,可以利用bedtools igv命令。

bedtools igv [OPTIONS] -i <bed/gff/vcf> #igv参数即可，[options]为图片保存路径

结尾介绍两个有瑞士军刀美誉的序列处理轮子seqtk （https://github.com/lh3/seqtk）和bedtools（http://bedtools.readthedocs.io/en/latest/）；小编经常使用，序列处理功能非常强大，这里简单介绍一下下~~

比如从fastq文件随机抽取部分reads

seqtk sample -s11 read1.fq 10000 > sub1.fq #-s随机数种子
seqtk sample -s11 read2.fq 10000 > sub2.fq #双端测序，需保持-s值一致

比如从bam中提取fastq

bedtools bamtofastq -i example.bam -fq example1.fq -fq2 example2.fq #生成双端reads,

当然这些工具功能的实现都可以自己编写脚本，如模仿bedtools bamtofastq命令

samtools view example.bam |perl -lane 'BEGIN{open IN1, ">example1.fq";open IN2,">example2.fq"}{if ($.%2==0){print IN1 "$F[0]\n$F[9]\n+\n$F[10]"}else {print IN2 "$F[0]\n$F[9]\n+\n$F[10]"}}'#利用perl 单行提取双端reads

值得注意的是bam文件一般按照染色体顺序进行排序，这两种方法都要求bam文件需要按照reads name进行排序。