手把手教你通過重測序尋找外源基因插入位點

在做轉基因動植物時，相信許多實驗室的小夥伴們都遇到過＂我轉成功了嗎？＂我轉到哪兒去了？＂這樣的困惑，那麼，你是否還在傳統的酶切，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輸出格式：

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,925snapshot 2.pnggoto Chr01:41,877,711-41,889,454Snapshot

然後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進行排序。