在做轉基因動植物時，相信許多實驗室的小夥伴們都遇到過＂我轉成功了嗎？＂我轉到哪兒去了？＂這樣的困惑，那麼，你是否還在傳統的酶切，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，925snapshot 2。pnggoto Chr01：41，877，711-41，889，454Snapshot
【左右滑動檢視完整資訊】
然後IGV工具欄Tools ——> Run Batch Script，就可以輕鬆批次擷取。如果基因組比較大，檢視區域也比較多，不太可能手動編輯如下指令碼，可以利用bedtools igv命令。
bedtools igv ［OPTIONS］ -i
【左右滑動檢視完整資訊】
結尾介紹兩個有瑞士軍刀美譽的序列處理輪子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進行排序。
好了，透過重測序尋找外源基因插入位點的方法就介紹完了，如果你還有不清楚的地方，可以給我們留言，小編看到會第一時間回覆你哦~~~
