转录组分析 | 使用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 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ，欢迎热爱写作的你一起参与！
腾讯云开发者
