早在本科學(xué)習(xí)比較基因組學(xué)期間就有意做一些知識(shí)分享设易，畢竟目前學(xué)得的七成生信知識(shí)是源自最早做知識(shí)分享的前輩們工三。計(jì)劃一直夭折，一方面是自覺實(shí)力有限炸宵，怕誤人子弟，一方面是有點(diǎn)空閑時(shí)間就想偷個(gè)懶谷扣。個(gè)人是希望做到干濕結(jié)合全面發(fā)展的土全，如果不想兩邊都學(xué)個(gè)半吊子，那必然得付出雙倍努力了会涎。我目前主要做倆物種裹匙，異源多倍體油菜和同源多倍體馬鈴薯，方向涉及多組學(xué)分析末秃、細(xì)胞遺傳和基因編輯概页。
回歸正題，基因編輯過的作物會(huì)有載體序列插入到植物基因組中练慕，確定T-DNA插入位點(diǎn)有重要用處惰匙，原理可以看一下這篇文章Illumina Sequencing Technology as a Method of Identifying T-DNA Insertion Loci in Activation-Tagged?Arabidopsis thaliana?Plants。下面介紹我是如何完成這項(xiàng)工作的铃将，有些內(nèi)容比如軟件安裝和參數(shù)設(shè)置项鬼，網(wǎng)上已經(jīng)有太多教程，這里就不啰嗦啦劲阎。
1绘盟、提取DNA二代測序，PE150悯仙，深度10X以上（太低可能檢測不到）龄毡。
2、過濾reads得到cleandata锡垄。
3沦零、創(chuàng)建樣本名文件samplename.txt，一個(gè)名字一行货岭。
4蠢终、運(yùn)行腳本
#!/bin/bash
#載體序列作為參考基因組，建索引
bwa index TDNA.fa
samtools faidx TDNA.fa
#讀入樣本茴她，寫個(gè)循環(huán)
cat samplename.txt | while read line
do
read1="${line}_1.clean.fq.gz"
read2="${line}_2.clean.fq.gz"
#bwa比對(duì)并samtools排序轉(zhuǎn)成bam文件
bwa mem -t 12 -R "@RG\tID:$line\tSM:$line\tLB:$line\tPL:ILLUMINA" TDNA.fa $read1 $read2 | samtools sort -@ 12 -o $line.sorted.bam
#samtools建索引，提取比對(duì)上的信息程奠，保存sam格式
samtools index -@ 12 $line.sorted.bam
samtools view $line.sorted.bam TDNA > $line.TDNA.sam
#提取比對(duì)上的reads的ID丈牢，根據(jù)ID從原始測序數(shù)據(jù)中提取這些reads
cut -f1 $line.TDNA.sam |sort|uniq > $line.TDNA.ID
seqtk subseq $read1 $line.TDNA.ID > ${line}_1.TDNA.fq
seqtk subseq $read2 $line.TDNA.ID > ${line}_2.TDNA.fq
#spades組裝到contig水平，contigs.fasta即為最終結(jié)果
spades.py --careful -1 ${line}_1.TDNA.fq -2 ${line}_2.TDNA.fq -o ${line}spades
done
5瞄沙、將contigs與載體序列blastn（圖1）己沛，比對(duì)不上的序列再和植物基因組blastn（圖2）慌核，即可找到插入位點(diǎn)。
精力有限申尼，難免出錯(cuò)垮卓，轉(zhuǎn)載請(qǐng)注明出處。有任何疑問师幕，歡迎交流討論粟按。
