转座子插入序列分析2-自制分析流程

我们先观察一下测序的结果，是否有一些什么规律，因为使用的靶向富集法的测序，我们使用了特殊序列将插入了转座子的部分钓了出来，然后进行的测序，所以理论上富集到的所有序列都应该存在一段与我们钓鱼序列互补的“靶点序列”。

我们可以看到，蓝色部分三个TTT前面的序列几乎相同，这是最典型的插入位点的特征，占所有插入位点特征的60%。同时也存在一些不典型的序列，没有TTT特征，事实上，依据公司给出的插入位点特征结果，仍存在40%的非典型插入位点结果（见下图），所以我们先对序列进行质控，使用cutadapt软件，将“钓鱼序列”视为接头序列进行删除。

我们想对序列5’端的进行质控，
我们先随机挑出来

cutadapt -g ATATCCCGCCAGGCCC -o ./004.newcatadaptQC/replicate1.1.fq Replicate1.fastq &

质控前

质控后

可以看见，带有典型“钓鱼序列”的序列已经经历了质控，而仍然有部分序列未正常质控，将这些序列拿出来blast，结果如下

我们先将这段序列质控掉，因为插入位点的因素，我们先尽量保留末尾的两个AT序列，如果影响了后续的BWA比对，再进行进一步质控。

这里面，大部分该序列已被删除，仍有少量序列未受影响，我们继续看这些序列有哪些特征。

我们仍然可以看到，这些序列的1-140bp相似度非常高，我们先删除完全相同的这些序列，后面根据比对结果再来调整，

比对

bwa mem -t 32 ../../human.ref.bwaindex/GCF_000001405.40_GRCh38.p14_genomic.fna 002.merged.3QC.fq > 002.1.merged.sam &

我们先去除所有没有比对上的序列。
提取文件的前六列

awk '{print $1, $2, $3, $4, $5, $6}' 002.1.merged.sam > 002.1.1-6merged.sam

去除表头

awk '!/@SQ/' 002.1.1-6merged.sam > 002.2.merged.sam

去除没有比对上的数据

awk '$4 != 0' 002.2.merged.sam > 002.3.merged.sam

然后，直接统计插入位点的数量，这样就能找到支持插入位点的reads数。

最后使用samtools+bedtools流程查看未找到位点的reads覆盖情况

samtools view -S -b merged.trimmed.QC.align.sam > merged.trimmed.QC.align.bam
bedtools coverage -a nomap.bed -b merged.trimmed.QC.align.bam > nomep.txt 

发现对应位点确实没有找到对应比对上的reads，接着逆推，提取位点下游15bp的序列，直接在fq文件中查找是否存在对应的reads。

bedtools getfasta -fi ../../../human.ref.bwaindex/GCF_000001405.40_GRCh38.p14_genomic.fna -bed nomap1-15bpseq.bed -fo nomap1-15bpseq.fa