专栏首页生物信息云转录组分析 | 使用trim-galore去除低质量的reads和adaptor

转录组分析 | 使用trim-galore去除低质量的reads和adaptor

我前面已经介绍了转录组分析中利用fastqc这个软件来查看测序质量【文章:转录组分析 | fastqc进行质控与结果解读】,通过分析结果报告,我测序的数据还是可以的,但不管怎样,还是要清除一些不好的reads。这里我用trim-galore去除低质量的reads和adaptor。

一.Trim Galore介绍

Trim Galore是对FastQC和Cutadapt的包装。适用于所有高通量测序,包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera 和smallRNA测序平台的双端和单端数据。主要功能包括两步: 首先去除低质量碱基,然后去除3' 末端的adapter, 如果没有指定具体的adapter,程序会自动检测前1million的序列,然后对比前12-13bp的序列是否符合以下类型的adapter: ● Illumina: AGATCGGAAGAGC ● Small RNA: TGGAATTCTCGG ● Nextera: CTGTCTCTTATA


当然,Trim Galore是可以自动检测adapter。

接下来我们看软件参数:

--quality:设定Phred quality score阈值,默认为20。我后面分析改成25,稍微严格一些。 --phred33::选择-phred33或者-phred64,表示测序平台使用的Phred quality score。具体怎么选择,看你用什么测序平台,这个在上一篇文章中的报告中就有【转录组分析 | fastqc进行质控与结果解读】。

具体是-phred33还是-phred64我在文章【生信中常见的数据文件格式】中有提到。

--adapter:输入adapter序列。也可以不输入,Trim Galore!会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个,也直接显式输入这三种平台,即--illumina、--nextera和--small_rna。 --stringency:设定可以忍受的前后adapter重叠的碱基数,默认为1(非常苛刻)。可以适度放宽,因为后一个adapter几乎不可能被测序仪读到。 --length:设定输出reads长度阈值,小于设定值会被抛弃。 --paired:对于双端测序结果,一对reads中,如果有一个被剔除,那么另一个会被同样抛弃,而不管是否达到标准。 --retain_unpaired:对于双端测序结果,一对reads中,如果一个read达到标准,但是对应的另一个要被抛弃,达到标准的read会被单独保存为一个文件。 --gzip和--dont_gzip:清洗后的数据zip打包或者不打包。 --output_dir:输入目录。需要提前建立目录,否则运行会报错。 -- trim-n : 移除read一端的reads

二.使用trim-galore去除低质量的reads和adaptor

首先,创建保存输出数据的文件夹。

 mkdir cleandata
 mkdir ./cleandata/trim_galoredata

然后接下来我们就可以处理数据了。下面命令是处理单个样本的案例,可以不先运行,只需要知道是什么意思。

trim_galore -q 25 --phred33 --stringency 3 --length 36  --paired CK-4_1.fq.gz CK-4_2.fq.gz --gzip -o ./cleandata/trim_galoredata/

参数解释前面已经介绍,这里提一下,我的测序是双端测序,我们前面12个文件对应6个样本,都是成对的。CK-4_1.fq.gz CK-4_2.fq.gz是一对,处理数据的时候就是一起处理。

但是我们有6个样本的数据,就需要写6个命令,好像有点不高效,所以我们可以写一个脚本,执行一次就行。

在当前路径下创建一个脚本文件trim_galore_batch.sh

vim trim_galore_batch.sh

输入下面内容:

#!/bin/bash
# This is for trimming a batch data
for i in CK-4 CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9
do
        trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired /data/RNAseq/${i}_1.fq.gz /data/RNAseq/${i}_2.fq.gz --gzip -o /data/RNAseq/cleandata/trim_galoredata
done

其实就是一个for循环:CK-4 CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9是文件名称的前部分,也就是样本名称,测序数据文件命名都是有规律的,认真看一下,就知道上面内容是什么意思。

保存脚本,然后执行脚本。

bash trim_galore_batch.sh

脚本没有问题的话,会看到类似下面的信息

(RNAseq1) [root@localhost RNAseq]# bash trim_galore_batch.sh
Multicore support not enabled. Proceeding with single-core trimming.
Path to Cutadapt set as: 'cutadapt' (default)
Cutadapt seems to be working fine (tested command 'cutadapt --version')
Cutadapt version: 1.18
single-core operation.
Output will be written into the directory: /data/RNAseq/cleandata/trim_galoredata/
AUTO-DETECTING ADAPTER TYPE
===========================
Attempting to auto-detect adapter type from the first 1 million sequences of the first file (>> /data/RNAseq/CK-4_1.fq.gz <<)

Found perfect matches for the following adapter sequences:
Adapter type    Count   Sequence        Sequences analysed      Percentage
Illumina        28671   AGATCGGAAGAGC   1000000 2.87
smallRNA        5       TGGAATTCTCGG    1000000 0.00
Nextera 4       CTGTCTCTTATA    1000000 0.00
Using Illumina adapter for trimming (count: 28671). Second best hit was smallRNA (count: 5)

Writing report to '/data/RNAseq/cleandata/trim_galoredata/CK-4_1.fq.gz_trimming_report.txt'

SUMMARISING RUN PARAMETERS
==========================
Input filename: /data/RNAseq/CK-4_1.fq.gz
Trimming mode: paired-end
Trim Galore version: 0.6.6
Cutadapt version: 1.18
Number of cores used for trimming: 1
Quality Phred score cutoff: 25
Quality encoding type selected: ASCII+33
Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 3 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 36 bp
Output file(s) will be GZIP compressed

Cutadapt seems to be reasonably up-to-date. Setting -j 1
Writing final adapter and quality trimmed output to CK-4_1_trimmed.fq.gz


  >>> Now performing quality (cutoff '-q 25') and adapter trimming in a single pass for the adapter sequence: 'AGATCGGAAGAGC' from file /data/RNAseq/CK-4_1.fq.gz <<<

这个时间很长很长。一个文件1个小时左右。慢慢等..................

运行结束后查看文件:

ll -h /data/RNAseq/cleandata/trim_galoredata

fq.gz格式文件是处理后得到的数据,如果还记得的话,前面我们的数据是27G,现在质控后只有22G的数据。txt格式文件是样品处理的结果报告,也包括软件运行的参数信息。下面是其中一个的结果。

SUMMARISING RUN PARAMETERS
==========================
Input filename: /data/RNAseq/CK-4_1.fq.gz
Trimming mode: paired-end
Trim Galore version: 0.6.6
Cutadapt version: 1.18
Number of cores used for trimming: 1
Quality Phred score cutoff: 25
Quality encoding type selected: ASCII+33
Using Illumina adapter for trimming (count: 28671). Second best hit was smallRNA (count: 5)
Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 3 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 36 bp
Output file will be GZIP compressed


This is cutadapt 1.18 with Python 2.7.15
Command line parameters: -j 1 -e 0.1 -q 25 -O 3 -a AGATCGGAAGAGC /data/RNAseq/CK-4_1.fq.gz
Processing reads on 1 core in single-end mode ...
Finished in 1496.55 s (62 us/read; 0.96 M reads/minute).

=== Summary ===

Total reads processed:              23,985,482
Reads with adapters:                 1,508,012 (6.3%)
Reads written (passing filters):    23,985,482 (100.0%)

Total basepairs processed: 3,597,822,300 bp
Quality-trimmed:               6,269,587 bp (0.2%)
Total written (filtered):  3,567,424,870 bp (99.2%)

=== Adapter 1 ===

Sequence: AGATCGGAAGAGC; Type: regular 3'; Length: 13; Trimmed: 1508012 times.

No. of allowed errors:
0-9 bp: 0; 10-13 bp: 1

Bases preceding removed adapters:
  A: 19.9%
  C: 34.7%
  G: 33.0%
  T: 12.4%
  none/other: 0.0%

Overview of removed sequences
length  count  expect  max.err  error counts
3  490193  374773.2  0  490193
4  129883  93693.3  0  129883
5  55548  23423.3  0  55548
6  30967  5855.8  0  30967
7  27693  1464.0  0  27693
8  27513  366.0  0  27513
9  25573  91.5  0  24770 803
10  24873  22.9  1  23690 1183
11  25827  5.7  1  24471 1356
12  25970  1.4  1  24553 1417
13  24143  0.4  1  23149 994
14  25530  0.4  1  24384 1146
15  24349  0.4  1  23344 1005
16  24383  0.4  1  23251 1132
17  24256  0.4  1  22978 1278
18  21920  0.4  1  21005 915
19  19694  0.4  1  18969 725
20  20484  0.4  1  19673 811
21  19751  0.4  1  18973 778
22  19049  0.4  1  18345 704
23  19704  0.4  1  18845 859
24  20170  0.4  1  19256 914
25  19510  0.4  1  18680 830
26  18750  0.4  1  18000 750
27  18256  0.4  1  17503 753
28  16883  0.4  1  16261 622
29  15923  0.4  1  15224 699
30  15198  0.4  1  14587 611
31  13894  0.4  1  13376 518
32  14121  0.4  1  13607 514
33  13453  0.4  1  12847 606
34  14484  0.4  1  13845 639
35  14362  0.4  1  13630 732
36  13907  0.4  1  13328 579
37  15268  0.4  1  14744 524
38  8918  0.4  1  8453 465
39  10640  0.4  1  10184 456
40  10771  0.4  1  10404 367
41  7766  0.4  1  7388 378
42  8181  0.4  1  7900 281
43  9762  0.4  1  9364 398
44  8479  0.4  1  8090 389
45  14727  0.4  1  14310 417
46  3156  0.4  1  2941 215
47  5179  0.4  1  4934 245
48  7148  0.4  1  6891 257
49  5975  0.4  1  5749 226
50  4647  0.4  1  4454 193
51  5267  0.4  1  5063 204
52  4253  0.4  1  4097 156
53  3958  0.4  1  3805 153
54  4345  0.4  1  4159 186
55  4843  0.4  1  4663 180
56  4152  0.4  1  3991 161
57  4613  0.4  1  4460 153
58  3645  0.4  1  3474 171
59  2935  0.4  1  2810 125
60  2315  0.4  1  2232 83
61  1453  0.4  1  1366 87
62  2271  0.4  1  2180 91
63  2382  0.4  1  2295 87
64  1702  0.4  1  1646 56
65  1394  0.4  1  1276 118
66  2765  0.4  1  2667 98
67  1689  0.4  1  1595 94
68  1880  0.4  1  1789 91
69  867  0.4  1  819 48
70  553  0.4  1  503 50
71  85  0.4  1  41 44
72  66  0.4  1  38 28
73  271  0.4  1  225 46
74  514  0.4  1  468 46
75  687  0.4  1  645 42
76  779  0.4  1  725 54
77  761  0.4  1  709 52
78  802  0.4  1  738 64
79  641  0.4  1  608 33
80  537  0.4  1  499 38
81  488  0.4  1  445 43
82  406  0.4  1  341 65
83  323  0.4  1  295 28
84  395  0.4  1  356 39
85  400  0.4  1  371 29
86  398  0.4  1  356 42
87  419  0.4  1  364 55
88  330  0.4  1  295 35
89  335  0.4  1  300 35
90  270  0.4  1  224 46
91  210  0.4  1  175 35
92  213  0.4  1  168 45
93  201  0.4  1  169 32
94  160  0.4  1  121 39
95  197  0.4  1  147 50
96  167  0.4  1  126 41
97  193  0.4  1  167 26
98  189  0.4  1  149 40
99  192  0.4  1  151 41
100  182  0.4  1  130 52
101  124  0.4  1  94 30
102  116  0.4  1  75 41
103  96  0.4  1  66 30
104  96  0.4  1  62 34
105  104  0.4  1  64 40
106  85  0.4  1  55 30
107  92  0.4  1  73 19
108  121  0.4  1  80 41
109  104  0.4  1  64 40
110  85  0.4  1  60 25
111  83  0.4  1  53 30
112  71  0.4  1  37 34
113  58  0.4  1  39 19
114  88  0.4  1  34 54
115  65  0.4  1  36 29
116  70  0.4  1  37 33
117  69  0.4  1  41 28
118  75  0.4  1  38 37
119  47  0.4  1  20 27
120  69  0.4  1  36 33
121  65  0.4  1  29 36
122  58  0.4  1  22 36
123  41  0.4  1  15 26
124  52  0.4  1  15 37
125  52  0.4  1  17 35
126  33  0.4  1  7 26
127  46  0.4  1  11 35
128  61  0.4  1  9 52
129  47  0.4  1  17 30
130  57  0.4  1  11 46
131  40  0.4  1  10 30
132  29  0.4  1  8 21
133  28  0.4  1  4 24
134  36  0.4  1  4 32
135  41  0.4  1  3 38
136  17  0.4  1  1 16
137  44  0.4  1  3 41
138  25  0.4  1  3 22
139  28  0.4  1  1 27
140  15  0.4  1  1 14
141  31  0.4  1  0 31
142  24  0.4  1  0 24
143  19  0.4  1  1 18
144  4  0.4  1  1 3
145  12  0.4  1  0 12
146  16  0.4  1  2 14
147  69  0.4  1  1 68
148  22  0.4  1  0 22
149  14  0.4  1  0 14
150  344  0.4  1  4 340


RUN STATISTICS FOR INPUT FILE: /data/RNAseq/CK-4_1.fq.gz
=============================================
23985482 sequences processed in total

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-09-16

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • TCGA数据库:生存分析

    本文介绍生存分析,其实,在R中,生存分析很简单,大家在网上能找到无数的文章。利用survival包就可以。就是按照下列公式就可以完成简单的生存分析。

    DoubleHelix
  • R语言基础绘图教程——第4章:面积图和饼图

    DoubleHelix
  • 细胞冻存实验

    实验开始前,将冻存管、 15 毫升离心管、移液管、移液枪、枪头等放入无菌超净工作台,以紫外线照射 30 分钟。采用通风机通风 3 分钟。以 75%酒精擦拭操作台...

    DoubleHelix
  • 将权重值映射到RGB 热力图

    最近做公司项目的时候遇到一个问题,怎样将1-20的权重值,用类似热力图的色值表示出来.起初直接rgb套用权重,出来的全都是灰度图.后来查了下资料,发现rgb三个...

    用户1220521
  • Visual Studio for Mac 2017 使用体验

    整个安装过程还是比较简单的,基本在官网安装包(约24MB),双击安装即可在线安装。

    Bug生活2048
  • 毕业设计项目,微博语料情感分析,文本分类

    微博的强大影响力已经深深的吸引了更多的人加入。而对微博的情感分析,不仅可以获取网民的此时的心情,对某个事件或事物的看法,还可以获取其潜在的商业价值,还能对社会的...

    机器学习AI算法工程
  • Linux CFS调度器之虚拟时钟vruntime与调度延迟--Linux进程的管理与调度(二十六)

    CFS为了实现公平,必须惩罚当前正在运行的进程,以使那些正在等待的进程下次被调度。

    233333
  • Oracle笔记之锁表和解锁

    开发过程经常遇到表被锁的情况,一般可能就是开发的修改数据库没提交事务,导致其他程序员不能再修改操作,这时可以用下面方法来解锁,这里主要涉及几张表

    SmileNicky
  • hibernate笔记(四)

    3) HQL查询, Hibernate Query language hibernate 提供的面向对象的查询语言。

    HUC思梦
  • TensorFlow 学前班

    本文我参加Udacity的深度学习基石课程的学习的第3周总结,主题是在学习 TensorFlow 之前,先自己做一个miniflow,通过本周的学习,对于Te...

    zhuanxu

扫码关注云+社区

领取腾讯云代金券