티스토리 뷰

앰플리콘 시퀀싱

Mothur 파이프라인 (pipeline)

메타지노믹스 메타지노믹스 2017.01.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




신고
댓글
  • 프로필사진 김효영 안녕하세요? 여기에 질문을 드리는게 맞는지 모르겠습니다^^

    Mothur 파이프라인중 silva.bacteria.fasta 파일이 필요한데요..
    혹시 다운로드 사이트 정보나 파일 제공이 가능하실까요?

    저에게 큰 도움이 될 것 같습니다. 정말 감사합니다~
    2017.01.23 17:24 신고
  • 프로필사진 메타지노믹스 메타지노믹스 Mothur 프로그램을 돌리려면 두개의 reference 파일이 필요한데요. 아래 링크에서 다운로드 받드실 수 있습니다.

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

    자세한 설명은 여기에 있습니다.
    https://www.mothur.org/wiki/Silva_reference_files

    2. RDP train set
    https://www.mothur.org/w/images/6/6c/Trainset14_032015.rdp.tgz

    설명은 여기: https://www.mothur.org/wiki/RDP_reference_files

    제가 데이터베이스 다운 받는 내용을 깜박했네요. 질문 주셔서 감사합니다. 위에 포스팅에 반영하겠습니다.
    2017.01.24 05:20 신고
  • 프로필사진 비밀댓글입니다 2017.01.24 14:47
  • 프로필사진 비밀댓글입니다 2017.01.24 22:36
  • 프로필사진 대학원생 안녕하세요.
    일루미나 Miseq Paired end data를 qiime으로 분석하다 여러 난관에 부딪혀 Mothur를 시도해보려하는데요.
    OTU pick이 아주 오래걸린다고 하셨는데 qiime closed reference OTU picking 과 비교했을 때 이보다 훨씬 오래 걸렸나요?
    교수님께서 주신 듀가 후덜덜해서 여쭤봅니다.

    미리 감사드립니다.
    그리고 포스팅들이 모두 따라하기 쉽게 잘 되어있어서 늘 감사드립니다.
    2017.03.03 14:10 신고
  • 프로필사진 메타지노믹스 메타지노믹스 필요한 시간은 데이터의 크기와 컴퓨터의 성능에 따라 크게 좌우 됩니다. 만약 데이터의 크기가 작다면 (1 기가 미만) 어떤 컴퓨터에서 작업 하셔도 오래 걸리지 않습니다. 다만 큰 데이터를 작업 하신다면 (4기가 이상) 저의 경험으로는 HPC(High performance computing)을 사용해도 1주 이상 걸렸습니다. 큰 데이터를 사용하신다면 QIIME이 좋은 방법이 될 수 있습니다. 예를 들어 제가 같은 크기 (4기가 이상)의 데이터를 Mothur와 QIIME에서 돌렸을 때, Mothur에서 1주일 이상 걸렸는데 QIIME에서는 1일 정도 걸렸습니다. 하지만 이건 저의 경험으로 말씀드리는 것이고, 가지고 계신 데이터의 Diversity 차이 등으로 인해서 저와는 다를 수도 있습니다. 도움이 되셨길 바랍니다. 2017.03.03 22:04 신고
  • 프로필사진 대학원생 감사합니다. 데이터량은 3기가이고 QIIME에서는 open reference에서 14시간, closed reference에서 5시간 정도 걸리더라구요. 꽤 오래 걸릴 것을 각오해야할 듯하네요.

    Illumina Miseq Paired end sequencing결과인데 barcode 시퀀스가 제거된 R1 fastq, R2 fastq file2개만 raw data로 받아서 mapping file이슈로 인해 split_libraries.py에서 더이상 진행이 되지않더라고요( barcode sequence부분 비운채로 mapping file을 만들고 validation은 -b option 으로 통과하더라도 split에서 mapping file error라고 계속 뜹니다.)

    그래서 mothur를 알아보고있었습니다. 감사합니다!
    2017.03.04 09:29 신고
  • 프로필사진 실례지만 질문 남기고 갑니다. metagenomics 공부 중에 있는 학부생입니다. 요즘 속해있는 연구실에서 mothur에 대해 공부하고 있는데 해당 글이 큰 도움이 되고 있습니다. 혹시 마지막 classify.otu쓰실 때 taxonomy parameter로 넣은 데이터가 뭔지 알고싶습니다. https://www.mothur.org/w/images/6/6c/Trainset14_032015.rdp.tgz 요 아이는 아닌거 같아서요..;; 2017.04.09 15:10 신고
  • 프로필사진 메타지노믹스 메타지노믹스 RDP trainset을 쓰는 단계를 깜박했네요. 위에 classify.seqs 단계에서 Trainset14_032015.rdp를 사용합니다. 마지막 classify.otu에서 사용된 파일은 전단계에서 생성된 파일입니다. 포스트 내용을 업데이트 하였습니다. 감사합니다. 2017.04.12 23:03 신고
  • 프로필사진 방황 정리 잘되어있는글 잘보았습니다 ㅎ
    많은 도움이 되었습니다.
    감사합니다

    몇가지 질문이 있는데 혹시 아시면 알려주시면 정말 감사하겠습니다.

    저희 실험실은 PGM으로 NGS 진행후에 CLC로 커뮤니티 분석을 진행하는데 이번에 네트워크 analysis를 위하여 Qiime 또는 mothur로 다시 분석중에 있습니다.
    근데 여기서 문제가 되는것이 PGM같은 경우 Barcode가 Trimming되어서 나오게 되는데, mothur를 구동하기 위해서는 oligo.txt를 만들어주어야 되더라구요 ㅠ 혹시 oligo.txt 없이 mothur구동할수 있는 방법이 있을까요??
    아시는점 있으면 알려주시면 정말 정말 감사하겠습니다 ㅜ

    일단 저는 리눅스를 하나도 모르는 상태에서 하나하나 따라가면서 진행중이긴 한데 ㅜ Miseq이나 454에 비해 PGM의 정보는 많이 부족하네여 ㅠ
    그렇다고 데이터를 Miseq으로 다시 돌릴수도 없고 그래도 메타지노믹스님의 친절한 포스팅에 조금이나마 힘이나네여 ㅜ
    2017.05.20 18:04 신고
  • 프로필사진 메타지노믹스 메타지노믹스 PGM 에서도 FASTQ로 파일을 받을 수 있지 않는지요? CLC가 있으시면 CLC에서 FASTQ로 변환하셔도 될것 같습니다. oligo.txt가 어떤 파일인지 제가 잘 모르겠네요. 혹시 조금 더 자세히 정보를 주실 수 있을까요? 2017.05.31 00:40 신고
댓글쓰기 폼
공지사항
Total
3,317
Today
8
Yesterday
20
링크
«   2017/06   »
        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  
글 보관함