转座子插入位点分析3------公司质控数据的自主分析

观察数据

数据修剪后的结果
可以看到，公司修剪后的结果均带有一个CGCCAGGCCC,而相应的，其后大多数带有“TTT”，这符合公司统计的插入位点碱基的统计结果，下图。因此，公司给出的结果是barcode+trimmed seq序列。

再看一下原始数据，找一下删除了那些序列
可以看见，修剪了CGCCAGGCCC之前的序列，同时修剪了所有不含CGCCAGGCCC的序列。在包含了CGCCAGGCCC的序列中，还修剪了3’端的序列，我们将这些序列比对一下，结果如下
可见这部分修剪掉的序列可能也是接头序列，被修剪掉了，而在我的之前的自主分析中，这部分序列被保留了，在BWT的比对算法中，如果匹配到正链中，这部分序列并不影响比对的起始位置，如果是匹配到负链，这部分序列会使其丢失匹配，这也是为什么我的数据中比对到负链的结果重复性很差的原因。

合并

由于这三个数据是重复测序数据，最终统计的插入位点是这三个数据的总和，所以直接将三个数据使用cat工具进行合并。

cat replicate1.trimmed.fq replicate2.trimmed.fq replicate3.trimmed.fq > merged.trimmed.fq

质控

由于公司给的数据中仍存在部分接头序列“CGCCAGGCCC”，这部分需要去除，否则会影响后续的比对结果
我们使用cutadapt对数据进行修剪

conda activate cutadapt
cutadapt -g CGCCAGGCCC -o merged.trimmed.QC.fq merged.trimmed.fq

运行日志

比对

使用bwa软件对将数据比对到人的参考基因组中。

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

比对结果的预处理

去除表头

awk '!/@SQ/' merged.trimmed.QC.align.sam > merged.trimmed.QC.align.1.sam

提取文件的前六列

awk '{print $1, $2, $3, $4, $5, $6}' merged.trimmed.QC.align.1.sam > merged.trimmed.QC.align.2.sam

去除没有比对上的序列

awk '$4 != 0' merged.trimmed.QC.align.2.sam > merged.trimmed.QC.align.3.sam

所有结果

最后，使用excel对MAP上的位置进行统计即可。