【生信学习笔记】fq数据过滤--接头问题

我的飞鸿踏雪泥

已于 2023-11-14 14:35:09 修改

阅读量556

点赞数

分类专栏：生信笔记文章标签：笔记

于 2023-11-14 14:18:41 首次发布

本文链接：https://blog.csdn.net/parkzx/article/details/134396919

版权

生信笔记专栏收录该内容

2 篇文章 1 订阅

订阅专栏

文章讨论了在处理fq测序数据时遇到的问题，发现某些Cleandata中的Readsid后面的index序列未去除。原因可能是插入片段过短导致测序通透。作者解释了文库结构和理想测序结果，指出实际数据中由于DNAinsert过短，可能导致测序数据结构异常，需要使用工具如cutadapt进行接头过滤。

摘要由CSDN通过智能技术生成

今天遇到一个客户提供的cut&Tag的fq测序数据，某公司测序提供的“Clean data”，好奇打开看了下，发现reads id 后面的index序列居然没有去除，如下图：

id部分保留index序列信息是不影响后续分析，但是第二行的序列尾部肉眼可见的index序列跟在后面，理论上来说，i5、i7index是用来拆分数据的，数据已经拆分成R1和R2了，不应该存在index序列在R1尾部。

通过对index的模式匹配发现，发现这份Clean data确实不够Clean。

那么出现这种情况出现的原因是什么呢？

从测序原理来看，应该是因为插入片段太短，导致测序测通导致的。

要想理解清楚这个问题，首先来了解下一般文库的结构（illumina）：

单index的SE测序

双index的SE测序

单index的PE测序

双index的PE测序

一般来说文库后的测序步骤应该是：

结合Read1 primer→开始测DNA insert→结合i7 index primer→开始测i7 index →结合i5 index primer→开始测i5 index→结合Reads primer→开始测DNA insert

（单双index的区别基本就是：测i5 index）

那么理想情况下我们得到的测序结果就是：

Read1+i7+i5+Read2

根据i7和i5的序列情况，我们就可以得到R1和R2的序列。比如PE150的数据，i7和i5一般是8bp的碱基，双index PE150测序的结果就是得到一整条长度为150+8+8+150的序列，根据这个结构拆分出R1和R2即可。

ps：i7和i5的另一个作用就是可以用来区分样本。

理想情况下，拆分出来的R1和R2应该就是DNA插入片段的两端的序列信息。

但是！！！

理想毕竟是理想！

压力给到下游的信息分析。

一般下游的信息分析进入分析的第一步就是数据过滤。

过滤的内容一般是过滤测序接头，低质量，含N多的Reads。

这里主要讲一下接头这个问题：

一般文库试剂盒会给出如下信息：

这里有两种引物的序列，我一开始的理解，这个接头应该指的是“测序引物”，因为理想情况下，DNA insert 在足够长的情况下，分析排除的应该是测序开始位置的错误，因为拍照、碱基识别算法、引物碱基结合等原因（我猜）（暗自揣测，为什么要给出多余信息）。

最近同事提出了这个接头的问题，深入思考了一下，只能说我当时太年轻，果然想法过于理想了。

从实际数据上来看这个DNA insert可能并没有我们理想的长，因为这种“随机打断”的方式很多，而且不同方式好像也有很大差别，长度好像并没有我想象中那么容易控制，原谅我不懂实验，所以不具体说明了。

那么从数据上来看，DNA insert太短就会导致将文库测通，即最终呈现的测序数据（生信分析拿到手的fq数据）结构就是：

R1：可能存在的Read1引物序列+真实的Read1序列+i7的index引物序列+i7

R2：可能存在的Read2引物序列+真实的Read2序列+i5的index引物序列+i5

那么真接头就要过滤掉所有的引物序列（暗自感慨:果然没有给多的信息）

ps：个人推荐用cutadapt进行过滤

使用网址：https://cutadapt.readthedocs.io/en/stable/index.html

我个人使用的参数如下：

#cutadapt -a 【i7 index primer】 -A 【i5 index primer】 -g 【Read1 测序引物】 -G 【Read2 测序引物】 -n 10 --max-n 0.01 -q 20,20 -m 15:15 -e 0.1 --report full -j 14 -o Clean_sample_R1.fq.gz -p Clean_sample_R2.fq.gz Raw_sample_R1.fq.gz Raw_sample_R2.fq.gz

#Nextera 文库
#Read1 primer：TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
#i7 index primer：CTGTCTCTTATACACATCTCCGAGCCCACGAGAC
#Read2 primer：GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
#i5 index primer：CTGTCTCTTATACACATCTGACGCTGCCGACGA

cutadapt -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA -g TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG -G GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG -n 10 --max-n 0.01 -q 20,20 -m 15:15 -e 0.1 --report full -j 14 -o Clean_sample_R1.fq.gz -p Clean_sample_R2.fq.gz Raw_sample_R1.fq.gz Raw_sample_R2.fq.gz

#具体参数的描述可以去官网查看：

我主要使用的参数说明如下：

cutadapt version 4.4

-a #Sequence of an adapter ligated to the 3' end (paired data: of the first read).
-g #Sequence of an adapter ligated to the 5' end (paired data: of the first read).
-A #3' adapter to be removed from R2
-G #5' adapter to be removed from R2
-n #Remove up to COUNT adapters from each read.
-e #Maximum allowed error rate (if 0 <= E < 1), or absolute number of errors for full-length adapter match (if E is an integer >= 1). Error rate = no. of errors divided by length of matching region. Default: 0.1 (10%)
--max-n #Discard reads with more than COUNT 'N' bases.
-q #[5'CUTOFF,]3'CUTOFF, --quality-cutoff [5'CUTOFF,]3'CUTOFF
-m #LEN[:LEN2] Discard reads shorter than LEN. Default: 0
--report #Which type of report to print
-j #Number of CPU cores to use.
-o #Write trimmed reads to FILE. FASTQ or FASTA format is chosen depending on input.
-p #Write R2 to FILE.