티스토리 뷰

앰플리콘 시퀀싱

Mothur 파이프라인 (pipeline)

메타지노믹스 2017. 1. 19. 00:33

Mothur 사용하기.


16s 앰플리콘 데이터를 분석하는 방법중 많이 사용하는 프로그램 중 하나인 Mothur를 사용려면 어떻게 해야 할까. 


Mothur 설치

Mothur의 가장 큰 장점이라면 '쉬운 설치'이다. 아래 웹사이트에 들어가서 자신의 컴퓨터의 운영체제에 맞는 파일을 다운로드 받아서 압축을 푼다. 설치가 끝났다!!!

https://github.com/mothur/mothur/releases

윈도우를 사용하는 경우라면 : Mothur.win_64.zip

맥을 사용하는 경우라며: Mothur.linux_64.zip

리눅스를 사용하는 경우라면: Mothur.linux_64.zip

을 다운로드 하면 된다.


데이터베이스 다운로드

Mothur에서는 두개의 reference database를 사용한다. 아래 링크에서 다운로드 받아서 압축을 풀면 된다. 

https://www.mothur.org/w/images/1/15/Silva.seed_v123.tgz

https://www.mothur.org/w/images/6/6c/Trainset14_032015.rdp.tgz


첫번째 파일은 Silva 데이터 베이스를 SEED Alignment를 이용해 Align한 것이다. 자세한 설명은 여기. 두번째 파일은 RDP에서 만들어 놓은 Training set으로 taxonomy assignment 를 할 때 사용한다. 자세한 설명은 여기 . 맥이나 리눅스를 사용하는 경우 여기를 참고하여 빠르게 다운로드/압축풀기가 가능하다. 


Mothur 실행하기

다운로드 받은 파일을 실행하면 된다. 윈도우 사용자라면 시작 > 실행 > cmd 를 입력하고 엔터, 를 통해 커맨드 창을 열 수 있다. (혹시 더블 클릭도 되나요???) 맥이나 리눅스 사용자는 자신의 운영체제에 맞게 커맨드라인에서 실행하면 된다. 커맨드라인을 사용하는 방법은 여기 링크를 참고 하면 된다.


Mothur SOP

Mothur 사용법은 웹사이트에 잘 안내되어 있다. 요즘은 대부분 앰플리콘 시퀀싱을 Miseq에서 돌리므로,, Mothur에서 제공하는 SOP 링크를 사용하면 된다. 개발자가 공식적으로 안내하는 방법이므로 이 방법으로 분석하면 나중에 논문쓸때도 편하다. (이 과정은 왜 이렇게 했는데!! 라고 꼬치꼬치 캐묻는 리뷰어의 공격 ㅜㅜ) 그런데 악!!! 영어!!! 영어 울렁증이 있는 사람을 위해 간단히 설명을 덧붙이면 아래와 같다.


1. 시작하기

Miseq에서 데이터를 받아 이미 Demultiplex를 끝냈다면, 파일이 F3D0_S188_L001_R1_001.fastq, F3D0_S188_L001_R2_001.fastq 이런 식으로 되어있을 것이다. 이 파일들의 정보를 주는 파일(stability file)을 만들어야 한다. 만드는 방법은 여기에 설명되어 있다. 간략하게 말해, 샘플이름, R1 파일이름, R2 파일 이름을 'tap'을 이용해 적어주면 된다. 파일 이름은 어느것이든 상관 없으나 이 예시에서는 stability.files 로 만들었다. 파일 내용은 아래처럼 만든다

sample1    sample1_R1_001.fastq    sample1_R2_001.fastq

sample2    sample2_R1_001.fastq    sample2_R2_001.fastq


2. Stability 파일을 이용해 시퀀스 파일들을 불러들인다.

아차, 먼저 Mothur를 실행한다. 그러면 커맨드라인이 아래처럼 바뀐다.

mothur >

그러면 명령어를 입력할 준비가 된 것이다.

그러면 명령어를 아래와 같이 입력한다.


mothur > make.contigs(file=stability.files, processors=8)


file은 위에서 만든 stability 파일을 넣어주면 된다. processors는 자신의 컴퓨터의 '코어 숫자'에 맞추어 넣어주면 된다. 예를들어 내 컴퓨터가 '쿼드코어'라면 4를 '듀얼코어'라면 2를 넣는다.


많은 경우에 아래와 같은 경고가 나온다. 그냥 무시한다. ㅋㅋ

[WARNING]: your sequence names contained ':'.  I changed them to '_' to avoid problems in your downstream analysis.


이렇게 하면 stability.trim.contigs.fasta 이 파일이 만들어 진다.


3. 시퀀스 길이로 필터링 하기

이 예제의 SOP는 V4 구역을 시퀀싱 한 것이다. 여러분이 다른 구역 (예를 들어 V1-V3)을 시퀀싱 했다면 그것에 맞추어 아래 숫자들을 바꾸어 주어야 한다. 길이가 275 bp보다 짧은 것은 버린다.


mothur > screen.seqs(fasta=stability.trim.contigs.fasta, group=stability.contigs.groups, maxambig=0, maxlength=275)


4. 반복되는 시퀀스를 버린다.

앗!.. 왜 반복되는 시퀀스를 버리는 걸까. 그렇게 되면 내가 원하는 시퀀스가 몇개 들어있는지 개수를 셀때 문제가 되지 않을까? 걱정마시라. 이 과정은 OTU 테이블을 만드는 과정에서 반복 시퀀스를 줄여 계산해야하는 경우의 수를 줄이는 작업이다. OTU를 결정하고 난 후에 원래의 개수를 다시 센다. 자 반복되는 시퀀스를 버리고 유니크 한 것만 모은다.


mothur > unique.seqs(fasta=stability.trim.contigs.good.fasta)


5. 카운트 테이블 만들기

카운트 테이블이 뭐지?? 카운트 테이블을 짧은 DNA서열들이 전체 중에 몇개씩 들어 있는지 세어 놓는 것이다. 무슨말이야!?!?! 이해가 어렵다면 그냥 컴퓨터의 계산을 빠르고 효율적으로 하기 위한 작업이라고 생각하면 된다.


mothur > count.seqs(name=stability.trim.contigs.good.names, group=stability.contigs.good.groups)


6. Align을 위한 준비단계

여기서 잠깐,, 우리는 레퍼런스(이미 알고 있는 시퀀스)에 우리의 샘플들을 정렬시킬 것이다. 그러기 위해서는 우선 레퍼런스 파일을 다운로드 받아야 한다. 그리고 우리는 V4구역만 필요하므로 V4에 해당하는 부분만 잘라낸다. 잘라내는 방법은? 컴퓨터로 PCR을 한다고 생각하면 된다. 여기서 나오는 숫자는 V4 구역이므로 여러분이 다른 구역을 사용했다면 이 숫자들이 바뀌어야 한다. Mothur SOP에는 연습용으로 만든 silva.bacteria.fasta를 사용했는데, 여기서는 실제 연구에서 사용할 수 있는 silva.seed_v123.align을 이용해 보자.

mothur > pcr.seqs(fasta=silva.seed_v123.align, start=11894, end=25319, keepdots=F, processors=8)


그리고 파일 이름을 바꾸어 준다.

mothur > system(mv silva.seed_v123.pcr.align silva.v4.fasta)


7. Align하기

자 이제 우리의 시퀀스를 레퍼런스에 align하자

mothur > align.seqs(fasta=stability.trim.contigs.good.unique.fasta, reference=silva.v4.fasta)


8. Summary file 만들기

mothur > summary.seqs(fasta=stability.trim.contigs.good.unique.align, count=stability.trim.contigs.good.count_table)


9. Align이 안된 녀석들 제거하기

Align이 안된 녀석들은 아마도 잘못된 녀석들일 것이다. 제거한다.

mothur > screen.seqs(fasta=stability.trim.contigs.good.unique.align, count=stability.trim.contigs.good.count_table, summary=stability.trim.contigs.good.unique.summary, start=1968, end=11550, maxhomop=8)


10. Overhang을 제거한다

여담으로... Mothur 개발자가 진행하는 워크샵에 TA로 간적이 있는데, 마지막 trump (트럼프)라는 부분에서 약간 당황해하며 (그때 트럼프가 대통령 후보였다) "이 (trump 라는) 명령어를 10년 전부터 사용했었다" 고 농담을 했었다. (지금 와서 보니 별로 웃기지 않은데 그때 참석자들이 모두 뒤집어 졌었음 ㅋㅋ)

mothur > filter.seqs(fasta=stability.trim.contigs.good.unique.good.align, vertical=T, trump=.)


11. 다시한번 반복 제거

mothur > unique.seqs(fasta=stability.trim.contigs.good.unique.good.filter.fasta, count=stability.trim.contigs.good.good.count_table)


12. pre-clust

mothur > pre.cluster(fasta=stability.trim.contigs.good.unique.good.filter.unique.fasta, count=stability.trim.contigs.good.unique.good.filter.count_table, diffs=2)


13. 키메라 찾기

mothur > chimera.uchime(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.fasta, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.count_table, dereplicate=t)


14. 찾은 키메라 제거

mothur > remove.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.fasta, accnos=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.accnos)


15. Taxonomy 찾기 

각 시퀀스가 어떤 녀석인지 찾아준다. 이 정보는 아래 단계에 이용된다. 

mothur > classify.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.count_table, reference=trainset14_032015.rdp/trainset14_032015.rdp.fasta, taxonomy=trainset14_032015.rdp/trainset14_032015.rdp.tax, cutoff=80)


16. 말이 되지 않는 녀석들 제거

16s 앰플리콘을 사용했으므로 진핵생물이 나올 수 없다. 만약 나왔다면 에러이므로 제거

mothur > remove.lineage(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.fasta, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.count_table, taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.rdp.wang.taxonomy, taxon=Chloroplast-Mitochondria-unknown-Archaea-Eukaryota)


17. [옵션] Mock 제거

만약 시퀀싱을 할 때 Mock community를 넣었다면 제거해준다. 만약 아니라면 이 과정은 skip 한다. 주의!! 만약 이 단계를 skip했다면 아래 파일 이름에서 pick.pick.pick 부분은 pick.pick으로 바꾸어야 한다.

mothur > remove.groups(count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.count_table, fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.fasta, taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.rdp.wang.pick.taxonomy, groups=Mock)


18. OTU 찾기

자!! 드디어!!! OTU 를 찾는다! 이 과정은 아~~~주 오래 걸린다. 만약 샘플이 작다면 [옵션1] 을 샘플이 크다면 좀 더 빠른 버전인 [옵션 2]를 사용한다.

[옵션1]

mothur > dist.seqs(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.fasta, cutoff=0.20)
mothur > cluster(column=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.dist, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.pick.count_table)


[옵션2]

mothur > cluster.split(fasta=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.fasta, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.pick.count_table, taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.rdp.wang.pick.pick.taxonomy, splitmethod=classify, taxlevel=4, cutoff=0.15)


19. 각 OTU에 얼마나 많은 시퀀스가 있는지 세기

mothur > make.shared(list=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.an.unique_list.list, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.pick.count_table, label=0.03)


20. Taxonomy 찾기

각 OTU가 어떤 녀석인지 Taxonomy를 찾아준다.

mothur > classify.otu(list=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.pick.an.unique_list.list, count=stability.trim.contigs.good.unique.good.filter.unique.precluster.denovo.uchime.pick.pick.pick.count_table, taxonomy=stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.rdp.wang.pick.pick.taxonomy, label=0.03)


자!!이제 끝났다!!! 와~ 엄청 길다!! Mothur의 단점은 모든 과정을 이렇게 직접 해야 한다는 것이다. 개발자도 이런 어려움을 느꼈는지 이 과정을 한번에 자동으로 끝낼 수 있는 방법을 마련해 놓았다.


자동화

자동화를 위해서는 위에 사용한 명령어들을 한개의 파일에 모아 stability.batch라는 이름으로 만든다. 실행은 아래와 같이 한다. (터미널에서 실행해야 함)

$ ./mothur stability.batch


위에서 설명한 파이프라인과 자동화를 위한 파일은 아래 링크에서 찾아 볼 수 있다.

https://github.com/metajinomics/mothur_pipeline




댓글
공지사항
최근에 올라온 글
최근에 달린 댓글
Total
Today
Yesterday
링크
«   2024/03   »
1 2
3 4 5 6 7 8 9
10 11 12 13 14 15 16
17 18 19 20 21 22 23
24 25 26 27 28 29 30
31
글 보관함