Raw data preprocessing

리눅스 (Ubuntu 16.0.4) 운영체제에서 분석한 내용을 토대로 작성한 글입니다.

(Illumina short read RNA-seq, Paired-end reads)

 

0. Raw data 

시퀀싱이 완료된 RNA-Seq raw data는 FASTQ 의 format 으로 되어있습니다.

 

주로 업체의 NGS report 에서 FASTQ file 을 받게 되고 주로 압축이 되어 있을 겁니다. 

 

ex) rawdata.fastq_1.tar.gz 또는 rawdata.fastq_1.gz

 

Paired-end sequencing data 는 샘플 하나당 forward 와 reverse 파일 한 쌍으로 제공됩니다. 

 

대부분 1과 2를 통해 구분하지만 업체마다 조금씩 다른 구분 방법이 있으니 미리 확인해 두는것이 좋습니다.

 

 

우선 압축부터 풉시다.

 

- 압축 해제하는 코드

 

.tar.gz 로 압축되어있는 경우

$ tar -xf rawdata.fastq_1.tar.gz

 

.gz 로 압축되어있는 경우

$ gzip -d rawdata.fastq_1.gz

 

 

압축을 해제 한후 rawdata.fastq 파일을 살펴보면, 눈썰미가 좋은 사람은 금새 4 가지 pattern 을 가진 line 의 반복인것을 알아 차릴겁니다.

 

FASTQ 파일에는 sequence 서열 뿐만 아니라 quality 정보도 포함하고 있는데

 

첫번째 줄은 “@” 로 시작하며 Sequence ID

두번째 줄은 Sequence 서열

세번째 줄은 “+” 하나만 있거나, 또는 그 뒤에 첫번째줄의 Sequence ID 의 반복

네번째 줄은 각 서열의 quality 를 나타내는 기호 (ASCII code) 로 이루어져 있습니다.

 

그렇다면 이 raw data 를 바로 활용할까요?

 

아닙니다. 

실제로, 거의 모든 연구에서 raw data 를 바로 사용하기보다는 preprocessing (quality control) 하는 과정을 거친 후에 이용합니다.

 

1. Preprocessing (quality control)

RNA-Seq raw data quality control 하는 tool 은 많지만 그중 제가 사용하는 건 Trimmomatic [1] (약 12,000 회 인용) 이라는 java 기반의 프로그램입니다.

 

Trimmomatic 홈페이지에 제공 되어있는 command

 

$ java -jar trimmomatic-0.35.jar PE -phred33 input_forward.fq input_reverse.fq \

$ output_forward_paired.fq output_forward_unpaired.fq \

$ output_reverse_paired.fq output_reverse_unpaired.fq \

$ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 \

 

저 같은 경우에는 더 빨리 하기위해 thread 수를 지정할 수 있는 –threads 옵션을 줍니다 (아마도 default는 4).

 

$ java -jar trimmomatic-0.35.jar PE -phred33 –threads 64 input_forward.fq input_reverse.fq \

$ output_forward_paired.fq output_forward_unpaired.fq \

$ output_reverse_paired.fq output_reverse_unpaired.fq \

$ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 \

 

또한 마지막에 있는 MINLEN:36 은 read trimming 과정중에 36bp 미만의 read 는 버리라는 명령인데,

저는 sequencing data 가 101bp 인 경우에는 MINLEN:36, 151bp 인 경우에는 MINLEN:50 을 줍니다.

* 중요 *

ILLUMINACLIP 뒤에 있는 TruSeq3-PE.fa 파일은 adaptor 파일입니다.

Trimmomatic 프로그램 내에 adaptor 디렉토리에 있습니다.

그래서 자기 sequencing data 에 맞는 adaptor 파일 경로를 지정해주면 됩니다.

Ex) ILLUMINACLIP:/home/lsg/toolbox/Trimmomatic/adaptor/TruSeq3-PE.fa

경로를 제대로 설정해주지 않으면 running 결과에 adptor 파일을 못 찾았다고 나옵니다.

~Results~

 

Trimmomatic 이 만약 잘 끝났다면, 아래와 같은 메세지가 출력될겁니다.

 

TrimmomaticPE: Started with arguments:

 -phred33 input_forward.fq.gz input_reverse.fq.gz output_forward_paired.fq.gz output_forward_unpaired.fq.gz output_reverse_paired.fq.gz output_reverse_unpaired.fq.gz

 

Using PrefixPair: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT' and 

'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'

Using Long Clipping Sequence: 'AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA'

Using Long Clipping Sequence: 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC'

Using Long Clipping Sequence: 'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'

Using Long Clipping Sequence: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT'

ILLUMINACLIP: Using 1 prefix pairs, 4 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences

Input Read Pairs: 23983840 Both Surviving: 18844813 (78.57%) Forward Only Surviving: 5046314 (21.04%) Reverse Only Surviving: 59830 (0.25%) Dropped: 32883 (0.14%)

TrimmomaticPE: Completed successfully

 

 

Both surviving 에 해당되는 (quality control 이 완료된) 결과는 output_forward_paired.fq, output_reverse_paired.fq 에 저장됩니다.

 

이 결과를 가지고 이후 분석 과정에 활용하면 됩니다.

  

------------------------------------------------------------------------------------------------------------

 

References

 

[1] Bolger, A. M., Lohse, M., & Usadel, B. (2014). Trimmomatic: A flexible trimmer for Illumina Sequence Data. Bioinformatics, btu170.