早在本科学习比较基因组学期间就有意做一些知识分享，毕竟目前学得的七成生信知识是源自最早做知识分享的前辈们。计划一直夭折，一方面是自觉实力有限，怕误人子弟，一方面是有点空闲时间就想偷个懒。个人是希望做到干湿结合全面发展的，如果不想两边都学个半吊子，那必然得付出双倍努力了。我目前主要做俩物种，异源多倍体油菜和同源多倍体马铃薯，方向涉及多组学分析、细胞遗传和基因编辑。
回归正题，基因编辑过的作物会有载体序列插入到植物基因组中，确定T-DNA插入位点有重要用处，原理可以看一下这篇文章Illumina Sequencing Technology as a Method of Identifying T-DNA Insertion Loci in Activation-Tagged Arabidopsis thaliana Plants。下面介绍我是如何完成这项工作的，有些内容比如软件安装和参数设置，网上已经有太多教程，这里就不啰嗦啦。
1、提取DNA二代测序，PE150，深度10X以上（太低可能检测不到）。
2、过滤reads得到cleandata。
3、创建样本名文件samplename.txt，一个名字一行。
4、运行脚本
#!/bin/bash
#载体序列作为参考基因组，建索引
bwa index TDNA.fa
samtools faidx TDNA.fa
#读入样本，写个循环
cat samplename.txt | while read line
do
read1="${line}_1.clean.fq.gz"
read2="${line}_2.clean.fq.gz"
#bwa比对并samtools排序转成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建索引，提取比对上的信息，保存sam格式
samtools index -@ 12 $line.sorted.bam
samtools view $line.sorted.bam TDNA > $line.TDNA.sam
#提取比对上的reads的ID，根据ID从原始测序数据中提取这些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即为最终结果
spades.py --careful -1 ${line}_1.TDNA.fq -2 ${line}_2.TDNA.fq -o ${line}spades
done
5、将contigs与载体序列blastn（图1），比对不上的序列再和植物基因组blastn（图2），即可找到插入位点。
精力有限，难免出错，转载请注明出处。有任何疑问，欢迎交流讨论。
