Entrez를 이용한 fasta 파일 다운받기

간만에 Biopython에 포함되어 있는 Entrez 함수를 이용하여 assceesion넘버로 fasta파일 다운받기를 해봤습니다.

git: https://github.com/gwlee/study/blob/master/entrez_access2fasta.py

python entrez_access2fasta.py accessionid 하면 {accessionid}.fasta파일이 생성됩니다.

참 쉽죠?

Biopython만 잘 사용하셔도 갱장한 것들을 하실 수 있으시고

그런 의미에서 다음번에는 좀더 재미진 내용으로 찾아오도록 하겠습니다. :)

출처: @sana_twice.09

bed파일을 편하게 다뤄주는 bedtools

bedtools를 설명하는데 swiss-army knife 만한게 없긴하지요
swiss-army knife 같이 다목적으로 적재적소에 쓰이는 좋은 물건 되겠습니다.

최신판은 2.29.2 되겠습니다. 당연한 얘기지만 상위버전으로 업데이트된 sub 명령어들의 옵션들이 추가되기 때문에 하위 버전 사용시 최근 메뉴얼 보시면 작동 안합니다. ㅋ
(웬만치 하위버전을 써야지.... 근데 이러는게 종종 있습니다. 그렇더라구요)

bed 포맷을 핸들링할때 굳이 스크립트 짜지 마세요
어설프게 짜고나서 나중에 후회마시고
그냥 bedtools 사용법을 남겨두고 써먹으시는게 서로 편합니다.

genome관련해서 사용할 만한 sub 명령어는

sort, merge, intersect, window 정도만 생각하고 있었는데.. 역시 coverage가;; 아놔.. 이거 겨우겨우 만들었는데...

역시 coverage라는 sub 명령어가;;;

여러분 힘들게 만들지 마시고 있는거 잘 활용하시면 

센스있는 분석쟁이로 가는 첫걸음입니다. :)

출처: @sana_twice.09

Long Read Assembler 설치 작업 로그

오랜만에 작업 로그용 글입니다. :)

Long Read(aka Nanopore)를 위한 assembler의 설치에 대한 로그로... 모 그렇게 자주 사용 될일이 없을것 같지만.. 그래도..

root권한 또는 sudo권한이 없는 상황을 가정하고 설치하는게...
나중에 편합니다. root권한 있으면 편하지만 나같은 쩌리한테 호기롭게 root권한이나 sudo를 부여할 이유가 있겠습니까? 그냥 없으면 없는대로 사는법도 알고 있어야... :)


canu (https://github.com/marbl/canu/releases)

$ wget https://github.com/marbl/canu/releases/download/v2.0/canu-2.0.Linux-amd64.tar.xz

$ tar -xvf canu-2.0.Linux-amd64.tar.xz

또는

$ git clone https://github.com/marbl/canu.git

$ cd canu/src

$ make -j <number of threads>


wtdbg2 (https://github.com/ruanjue/wtdbg2)

$ git clone https://github.com/ruanjue/wtdbg2

$ cd wtdbg2 && make


Raven (https://github.com/lbcb-sci/raven)

$ git clone --recursive https://github.com/lbcb-sci/raven.git raven

$ cd raven && mkdir build && cd build

$cmake -DCMAKE_BUILD_TYPE=Release .. && make

$ ./bin/raven

단, raven은 cmakr 3.9이상이 필요합니다. cmake 설치는 아래에 따로..


Racon (https://github.com/lbcb-sci/racon)

$ git clone --recursive https://github.com/lbcb-sci/racon.git racon

$ cd racon

$ mkdir build

$ cd build

$ cmake -DCMAKE_BUILD_TYPE=Release ..

$ make

racon의 경우 raven이 아닌 miniasm_and_minipolish.sh 작업시 racon을 찾아 해매서 racon 설치도 진행하였습니다.


flye (https://github.com/fenderglass/Flye)

$ git clone https://github.com/fenderglass/Flye

$ cd Flye

$ python setup.py install --prefix=/path/to/install/

또는

$ python setup.py install --user

※ --user 라는 옵션이 갱장히 편합니다. 대신 나만 됩니다.





cmake (https://cmake.org/)

$ wget https://cmake.org/files/v3.10/cmake-3.10.3.tar.gz

$ /bootstrap --prefix=/path/to/install/

$ make

$ make install

※ prefix를 설정하지 않으면 /usr/bin 모 이런데에 설치 되므로 설치가 제대로 되지 않기 떄문에 prefix를 설정하는것이 정신건강에 이롭습니다. :)

출처: @sana_twice.09

Benchmarking of long-read assemblers for prokaryote whole genome sequencing

나노포어는 현존하는 시퀀서중에 가장 긴 서열을 뽑아내는 시퀀서임에는 그 누구도 부인하지 못할것입니다. 근데.. 생산된 리드의 각 base의 phred score를 보자면.. 왜 갑자기 눈에서 물이나오는 이유는 왜때문일까요?
(그렇지만 저는 de-novo할때 보수적인 그룹이 아니라면 나노포어를 권장하는건 비밀..)

여하튼.. 현재 나노포어 어셈블리 용으로 이런저런 어셈블러가 판치고 있는 난세에 누가누가 좋은지 확인하는 작업을 해서 투고하신분이 나타나셨습니다.
제목도 정직합니다. 단, prokaryote대상입니다.
Benchmarking of long-read assemblers for prokaryote whole genome sequencing

prokaryote에서도 개판이면 굳이 사용할 이유가 있겠느냐? 주의 되겠습니다.
일단 가장 좋은것은 모르겠지만 최악은 걸러내야 해야 시간 낭비, 전기 낭비 하지 않지 않겠습니까?

현재(aka 당시에) 돌려볼 수 있는 7개 어셈블러 (Canu, Flye, Miniasm/Minipolish, NECAT, Raven, Redbean, Shasta)의 성능을 비교 평가 했습니다.
어셈블리의 정확성은 당연하고, prokarypte다 보니 circularisation도 중요하고, 계산시 사용되는 리소스와 분석 시간등을 평가했다고 합니다.

아름다운 figure는 상단에 링크된 논문에서 감상하시면 되고,
canu는 그나마 볼만한 서열들을 제공해줬고
flye는 canu다음으로 괜찮은 서열로 어셈블리 했다고 합니다.
redbean(wtdbg2) 과 shasta는 계산 리소스와 분석 시간에서는 효율적이었지만 결과는 그다지 효율적이지 않았고 하네요.

그래서 종합해서 논문에서 결론을 냈는데
모.. de-novo aseeembly 해보신분이라면 알고계시다 싶이.. 다들 장단점이 있었고, 원탑인 어셈블러는 없었지만 그 중에서 Flye, Miniasm / Minipolish와 raven이지 않나 싶다고 하네요

Flye는 믿을만한 서열을 제공했고(low depth에서도 나름..)
Miniasm / Minipolish는 circularisation이 좋았고
raven은 identity가 낮은 read set들에서 tolerant가 있었다고 합니다.

역시 최적의 어셈블리를 위한 정도는 당신이 사용 가능한 리소스를 동원해서 다양하게 돌려보고 비교한 게 킹왕짱이지 남의말 믿고 쓰면 너만 바보 되고
개발자님들인 이런 상황이니 개발좀 굽신굽신 :)

출처: @sana_twice.09

DTC는 얼마나 벌고있을까?

SARS-CoV-2로 어수선한 요즘이지만 이전부터 유전체 검사중 DTC(Direct To Consumer, 소비자 직접 의뢰) 상품/서비스들은 과연 돈은 벌고 있는지 공시자료에 그 정보들이 있지 않을까해서 한번 찾아봤습니다.

확인 기준은 현재 2019년도 사업보고서가 올라오지 않아서 2019년 9월이 제가 확인 가능한 최신 자료입니다.

	대분류	상세분류	단위: 백만원
마크로젠	용역매출	DNA Sequencing Microarray	77,356 3,855
테라젠이텍스	유전체사업(용역)	제노맘(NIPS)	1,140
		제노브로(착상전배아검사)	351
		유전체분석(기타용역서비스)	102
		유전체분석(Hellogene/Total Omics)	9,496
		진스타일SKIN	187
디엔에이링크	개인유전체분석	DNAGPS	620
랩지노믹스	용역	일반진단	10,839
		분자진단(유전체)	8,024
녹십자	유전자분석	유전체분석용역등	8,547
	임상검사분석	임상시험검사분석서비스	813
EDGC	유전체진단서비스매출	비침습산전진단	1,925
		신생아유전질환스크리닝	462
		안과질환유전체분석	6
		진료과별맞춤(gene2me)	459
		DTC(gene2me)	98
		유전성유방/난소암진단(BRCEARE)	40
		피부진단검사(Skincare)	12
		생체나이진단검사	43

마크로젠도 DTC사업을 진행하고 있지만 공시자료로는 어떻게 얼마나 매출을 올리는지 알기가 어렵네요. 랩지와 녹십자도 좀 구분하기는 어려웠습니다. 테라젠, 디링과 EDGC는 그래도 잘 구분해 놓아서 DTC 사업으로 어느정도 매출을 올리고 있는지 대략적으로 알 수 있었습니다. :)

사실 DTC는 순수히 고객이 요청해서 받을 수 있는 DTC의 매출을 정확히 확인 할 수 있었던 기업은 그나마 구분을 잘해놓은 테라젠과 EDGC정도? 디링의 DNAGPS는 DTC와 의료기관에서 진행가능한 유전체 서비스가 섞여있는것으로 알고 있어서...

갑자기 나는 왜 DTC 사업을 영위하고 있는 업체들의 DTC 매출을 궁금해 했는가...?
(aka 아 저 인간이 또 무슨 헛소리를 할려고 이걸 조사했지? ㅋㅋ)

출처: @sana_twice.09

현재 SARS-CoV-2 Tree를 그리면 어떻게 나올까

2월 28일 기준 2월의 마지막 주말을 맞아 gisaid에 몇개나 업로드 됐는지 확인하러 들어 갔는데 200개가 넘었네??
(2020년 2월 29일 현재 234개임)

그래서 tree 그려보면 어떻게 나올까 궁금해서 한번 그려봤습니다.
근데 200여개 중에 full length로 보이는 거는 대략 164개 정도

그래서 164개 골라내고 SARS 4개 서열과 MERS 2개 서열을 함께 align하고 tree를 그냥 후딱 그려보았다.
(집에서 그냥 작업용으로 사용하는 PC가 i3에 그냥 서류작업에 샘플 테스트는 돌릴 수 있는 정도지 무엇인가 작업하기에는 어렵지 않겠습니까? 그래서 다음과 같이...
최소한의 부하가 걸리지 않는 옵션으로...)

Align: mafft in.fa > out.fa (정확도를 높이는 작업은 진행안했다. 집 컴퓨터로 하려니 안끝나서)
Phylogentic Tree: MEGA X (Maximum Likelihood, General Time Reversible model 사용 했고, 당연히 Bootstrap 100이도 해보려고 했는데 안끝나서 결국 컴터를 껐다는 Orz)

근데 서열중에 Pangolin이라는 단어에 2017이라는 숫자가 있어서 뭐지? 했는데
이게 바로 그 천산갑  ㄷㄷㄷ
천산갑을 벌써 시퀀싱해서 올렸나? 근데 샘플링날짜가 2017? 모지.. 여튼


여튼 그냥 집에서 대충 취미삼아 그려본거...
우리나라에서 올라온 서열들은 7개가 있는데 SNU01, KCDC12, KCDC03/05/06/07/24 3그룹으로 나눠지는 듯한.....

이거 가지고 알 수 있는것은 7명에서 얻은 샘플이 gisaid에 올라와 있고 다른 서열들과 비교해보니 이렇다더라 정도...

그냥 가볍게 보시면 될것 같습니다. :)

GISAID 164ea + Other

GISAID 164ea

ps. 음... 위의 이미지를 저도 한번 봐봤는데 잘 안보이네요..

그렇게 꼭 보고 싶으시면 메일 보내주시면 대단한것도 아니니 164개와 SARS/ MERS 서열 파일 보내드리겠습니다.

출처: @sana_twice.09

2019-nCoV Tree 그려보기 -End-

지난번 2번째 글에서는 NCBI에서 genbank파일을 다운로드 받아서 python script로 어쩌구 저쩌구하면 2019-nCoV 서열을 모을 수 있다고 했는데요...

2019-nCoV관련 더 많은 정보들을 확인하시려면 gisaid라는곳에서 서열 을 받으시면 되겠습니다.

gisaid.org

현재(2020년 2월 8일) 76개 서열이 업로드되어 있고 그중 complete genome이 아닌 몇개가 있어서 대략 70개의 서열이 업로드되어 있다고 보시면 되겠습니다.

(아.. 아래 화면은 당연히 회원가입 해야 확인 할 수 있습니다.)

여기서 몇몇 1번 서열들을 다운로드해서 테스트 해보도록 하겠습니다. (최근에 KCDC에서 발표한 서열도 포함하였습니다.

그럼 MEGA-X를 수행하기 전에 서열 정렬 프로그램으로 mafft를 사용하도록 하겠습니다.
(muscle도 좋은데 mafft가 더 빨라서... 빠르다고 좋은건 아니지만 어차피 dna서열 MSA방법이 손꼽는 특정 프로그램이 없는 관계로.. clustalw, clustalo, muscle, tcoffe, mafft 정확히는 protein 서열 정렬이 주특기이고 protein 서열 정렬이 더 의미가 많...)

여하튼.. 그래서

>mafft.bat --auto 2019-nCoV_10ea.fasta > 2019-nCoV_10ea.mafft.fasta

한 결과를 MEGA-X에서 열어보면 이렇게!!

몬가 잘 서열이 잘 정렬 된것처럼 보이죠

Phylogenetic Analsis를 클릭하고, Confirmation창이 하나 뜨는데 저는 No 선택합니다. 이 서열은 전체 genome이지 coding하는 서열이 아니라서.. (물론 coding되는 서열들입니다. ㅎㅎ)

그리고 난 후

MEGA-X 실행 화면에서 PHYLOGENY 선택하고,

5가지방법중 맘에 드시는거 선택해서 작업하시면

진행바가 죽죽 진행되면서

짠하고 다음과 같이 Tree가...

이렇게 그려집니다.

다음번에 기회되면 Tree 그릴때 사용되는 방법에 대해서 공부해서 작성해보는 기회가 있기를 바래보며... 

2020년 새해 벽두 (설날기준)부터 전세계를 공포로 몰아넣고 있는 2019-nCoV Tree그리기를 마무리 하도록 하겠습니다. :)

출처: @sana_twice.09

2019-nCoV Tree 그려보기 -2-

2019-nCoV Tree 그려보기 1편에서 어라?
실망하셨던분들을 위해 준비한 2편 되겠습니다.

일단 일전에 말씀드린것과 같이 NCBI 홈페이지가서 다음과 같이 genbank파일을 입맛에 맞게 받으시면 되겠습니다.

refseq말고 다 받고 싶으시면 다 받으셔도 되요.

[Create File] 버튼을 클릭하시면 sequences.gb? 라는 파일로 다운로드가 될것입니다.
(사실 fasta파일로 받으셔도 상관은 없는데... genbank파일로 받으시는게 나중에 더 많은 일을 하실 수가 있으십니다.)

(20년 2월 3일 기준 coronavirus nucleotide 서열 다받으니 대략 300M정도 나왔습니다.)

그럼 이 파일가지고 어쩌라고?
라떼는 fasta파일 가지고.. 어? 마? 그랬는데?
그래서 준비했습니다. gb파일을 fasta파일로 만들어주는 바로그 스크립트!

import os,sys
from Bio import SeqIO
inFile = 'sequences.gb'
for rec in SeqIO.parse(inFile, "genbank"):
seq_id = rec.id
seq_name = rec.name
seq_desc = rec.description
seq_seq = rec.seq
print '>{}|{}\n{}'.format(seq_name,seq_desc,seq_seq)

이렇게 하면 다운바은 sequences.gb파일의 nucleotide 전체 서열을 fasta파일로 만들어주게 됩니다. 만약 나는 다 필요없는데... 하시는 분은
print 전에 if문 넣어서 특정 seq_name의 서열들만 혹은 seq_desc안에 특정 문자열.. 예를 들어 wuhan과 같은 내용이 있는 서열들만 선별해서 fasta파일로 만들 수 있습니다.
참 쉽죠?

그럼 다음 이시간에는 MEGA-X를 사용하는 방법을 얘기해보도록 하겠습니다. 제발~

출처: @sana_twice.09

2019-nCoV Tree 그려보기 -1-

2020년 정월 초하루가 얼마 지나지 않은 나날 동안 한중일이 아닌 전세계가 코로나 바이러스 때문에 홍역을 치르고 있습니다.

초기 방역로 이렇게 급속도로 전세계로 퍼져나가 감염자 및 사망자가 증가하고 있는 가운데 이 짧은 기간 동안에 다양한 연구 결과들이 나오고 있는데...

신종코로나 바이러스 논문들을 보고 있노라면 정렬을 왜그리 해대는지 그리고 무슨 나무를 왜그렇게 그려대는 걸까? 라는 궁금증을 가지시는 분들이 있을 것이라고 생각해서 그런게 무엇인지 그리고 나는 그려볼 수 없는지 하시는 분들을 위해서 MEGA를 이용한 corona virus phylogenetic tree를 한번 그려보고자 합니다.

이거대로 따로하면 당신도 tree 전문가(는 아니고 그냥 tree draw skill +1)

연구자들이 NCBI에 서열들이 올라와 있다고 얘기를 하고 서열들을 서로 비교한다는데
그럼 NCBI는 어떻게 접속하고 서열들은 어떻게 받지?

NCBI는 구글에 물어봐도 되고 네이버에 물어봐도 됩니다.
https://www.ncbi.nlm.nih.gov/

여기에 그냥 검색하듯 coronavirus를 검색하면 다음과 같이 나오는데..

이걸 어떻게 받나?
하나하나 집념으로 다운받으셔도 되고요...

다음과 같은 방법도 있습니다.

아까의 화면에서 Viruses, genomic, RefSeq 항목을 클릭해주면 해당 항목을 만족하는 서열들만 선택이 됩니다.

그리고 오른쪽 상단에 있는 Send to: 라는 항목을 클릭하면...!!
다음 이미지에 나와 있는것처럼 다운 받으시면 되겠습니다.

다운로드한 genbank파일을 parsing해서 비교해보고 싶은 부위의 서열을 득해서
MEGA라는 프로그램을 다음 url(https://www.megasoftware.net/)에서 다운로드 받아 설치한 후 입력한후 서열 정렬 및 Tree 작업을 해주면...

요렇게 이뿌게 Tree를 만들어 줍니다.

이렇게 포스팅을 끝내면 저게 모야? 라고 쌍욕을 날리실겁니다!

그래서 Tree 그릴때 필요한 서열을 수집하는 작업과 관련 스크립트에 대해서도 자세히 적어보는 글은 조만간 투척하도록 하겠습니다. :)

출처: @sana_twice.09

genpept 파일을 parsing하고 싶을때

genpept는 무슨 파일인고?

자세한 설명은 >여기<

간단히 얘기하자면 예전에 사용하셨던 genbank 파일 포맷입니다.
이름만 바꾼건지 모 그렇습니다. ㅎㅎ
biopython에서도 genpept가 아니라 genbank를 사용해서 접근하면 잘 parsing됩니다.

import os,sys
from Bio import SeqIO
try:
    inFile = sys.argv[1]
except:
    print ''
    exit(1) 

for seq_record in SeqIO.parse(inFile, "genbank"):
    print (seq_record.id)
    print (seq_record.name)
    seq_anno = seq_record.annotations
    print (seq_anno['accessions'])
    print (seq_record.seq)

이렇게 코드 짜서 사용하시면 됩니다. annotations안의 정보는 dict형식이라서 dict사용하는방식으로 확인 할 수 있습니다.

간만에 gp파일 다운받아서 정리하는 김에 biopython의 genbank 사용법 정리해보았습니다.  :)

크.. WSL이 있으니 참 편하긴 합니다.

출처: @sana_twice.09