너도 AB1파일가지고 pdf 만들수 있어 (1)

국내 회사에서 일반 시퀀싱 즉 Sanger Sequencing을 맡기시면

시퀀싱되어 바이너리 정보로 저장되어 있는 ab1파일, 서열정보가 들어있는 fasta파일과 함께 Applied Biosystems 로고가 찍혀있는 pdf파일 혹은 시퀀싱을 요청한 회사의 로고가 찍혀있는 pdf파일을 받아 보실 수 있으실겁니다. 

자세히 보다보시면 서열과 plot의 위치가 맞는거 같으면서도 안맞는거 같기도하고..

대놓고 틀린건 아닌데 그렇다고 안대놓고 맞는건 아닌거같고..

그러면 이 파일은 대체 어떻게 만들어지는걸까?

그럼 나는 저런거 못만드나?

그럴리가 있겠습니까?

그래서 이 글의 제목이 "너도 AB1파일로 pdf 만들 수 있어" 아니겠습니까?

예전에는 ab1.py? ab1_parse.py? 어떤 용자께서 제작하신 3rd party파일을 python에 import하여 사용하였었는데

이번에는 biopython에서 제공되는 ab1파일을 핸들링 하는 기능을 사용해서 작업해보려고 합니다.

ab1파일을 parsing하려면 ab1파일이 어떻게 생겼는지 알아야 하겠죠?

그래서 ab1파일에 어떤 항목에 어떤 내용들이 담겨져 있는지 확인을 해보시기 바랍니다.

>ABIF_File_Format.pdf<

다음번에는 ab1파일안에 들어있는 내용들을 확인해보는 시간을 가지도록 하겠습니다.

진짜로~

제발~

출처: @sana_twice.09

WGS는 과연 의료비용을 낮출수 있을까

Rapid whole-genome sequencing decreases infant morbidity and cost of hospitalization

본 논문은 WGS이 의료 비용 절감과 환우의 치료에 대해 긍정적인 결과에 대한 내용을 담고 있는데 솔까말 시퀀싱 가격이 NovaSeq이 나오면서 값싸졌다고는 하지만 의료 현장에서 WGS가격이 수십만원대로 서비스할 수 있는 것은 아니지만...

환자에게 들어가는 의료 비용을 생애 전주기에 걸쳐있다는 관점으로 바라보면, 초기에 WGS를 수행해서 병의 정확한 이유를 가려내어서 알맞은 치료나 관리를 해주면 초기에 다소 비싸보이는 WGS를 하는것이 WGS를 하지 않아 적시에 필요한 치료과 관리를 하지 못했을 때 보다 전체적인 의료비용도 낮추고 환우들이나 그의 가족들에게도 물심양면적으로 긍정적인 효과가 된다는 내용 되겠습니다.

San Diego의 Rady Children 's Hospital에서 진행되고 있는 연구로 WGS데이터를 바탕으로 질병 예측과 이유의 정확도를 높여 더 나은 치료 결과를 기대 할 수 있는지를 보고자 하는 임상실험 결과 중 초기 단계의 결과입니다. 

무슨 임상실험을 2050년까지 한다는건지.. 그러나 그정도로 데이터를 쌓아놓고 확인해봐야 NGS를 의료계에서 인정하겠다라는 철저한 준비를 하겠다는 의지라고 볼 수도 있겠지만.. 30년동안 데이터 모으는 동안 또다른 혁신적인 실험기법/방법이 나오면 도루묵일수도... 여튼.. 

그래도 의료에 직접적으로 영향을 끼칠 수 있는 방법이다보니.. 

그리고 이러쿵 저러쿵 하더라도 중요한 점은 FDA 승인하에 의뤄지는 WGS기반의 진단 방법 실험 결과이고 그 결과가 "가치 있음" 이라는데에 의의가 있겠습니다.

2016년/2017년 이유를 알지 못하는 유전 질환을 앓고 있는 영아 48명 중 유전자 검사에 대해서 동의한 42명의 영아들에 대해서 (필요시 영아와 함께 부모또는 형제자매들도 검사) 유전자 검사를 진행하였습니다.

rWGS의 임상적 효과를 확인하기위해 standard test인 염색체 마이크로어레이(CMA)와 함께 진행하여 비교 한 결과, CMA방법으로는 42명 환자중 4명에 대해서 병인에 대해서 진단 할 수 있었으나, rWGS의 경우 18명(?19명?)에 대해서 병인 진단을 받았고 위양진단은 없었음을 확인하였다고 합니다. 

그리고 유전질 환을 진단받은 18명의 영아중 13명에 대해서 소위 이야기하는 정밀/ 맞춤의학이 도입되었을 때 볼수 있는 효과는 Table 2,3에 설명되어 있습니다. 물론 재현성이 있는 반복실험 결과도 없고 코호트내 인종이 caucasian, Hispanic/Latino로 편중되어 있고, n수 도 수십명으로 적어서 rWGS를 적극적으로 의료보험으로 적용시키기에는 제한적이긴 하나 rWGS가 현재 적용되고 있는 기준 방법인 CMA방법보다는 더 나은 선택지를 제공할 수 있고 그 효과는 무시할 수 없다는 근거가 이 논문을 시작으로 속속 뒷받침 해줄 것이고 원인을 명확하게 판단 할 수 없는 영유아의 유전질환 진단에서 rWGS가 조만간 단계적으로 도입되는 날이 얼마 남지 않은것 같아보입니다.

사실 이 논문의 배경은 미국이라 우리나라에는 적용하기 어려운 부분이 많아보이는데 국민건강보험기금이 화수분도 아니고 이제는 조금더 효율적으로 사용하려면 신기술에 대해서 적극적으로 검토하고 도입을 주도적으로 진행해야 적은 기금으로 효과를 극대화하고 필요한 많은 이들에게 보장을 해 줄 수 있지 않을까하는 생각을 들게하는 연구였습니다.

출처: @sana_twice.09

Uterine Microbiota의 이해

읽는다 읽는다 하면서 이제서야 정리하는 내용...

Uterine Microbiota: Residents, Tourists, or invaders?

논문좌표 : https://www.frontiersin.org/articles/10.3389/fimmu.2018.00208/full

microbiota라고 하면 대부문 gut microbiota를 떠올리지만 인체에는 다양한 외부 기관들이 있는바 거기에 각각 독특한 microbiota를 가지고 있다는 사실!!

그중에서도 오늘은 글쓴이가 한때 업무를 보았던 NIPS와 밀접한 자궁내 미생물총에 대해서 리뷰된 논문이 있어서 한번 읽어 보았습니다.

설마 재미가 있을리는.... 

업무와 관련있는 내용을 논문이 나온지 2년이나 지난 지금에서야 읽고 정리하는 역설적인 상황!! 

어쨌거나 논문을 들여다보자면...

일단 리뷰지인만큼 미생물총, 그중에 자궁 미생물 총에 대해서 배경을 설명해주고 있습니다. 1세기 전만해도 자궁은 무균지대라는 가설이 지배적이었다고... 질은 외부세균으로부터 자궁을 보호하는 역할을..  (아니 무슨 말을 하는건지.. 밁칡짉안곩섥는... 사실 먼 훗날 현재는 사실로 받아들이는것들이 헛소리라고 밝혀지는것도 있을테니....)

그래도 최근에는 culture방법을 통해서 저 말도안되는 가설을 밀어내고 있고 현재는 NGS로 다양한 연구를 통해 자궁내 미생물총에 대한 기능에 대해서 알아가고 있습니다.

다양한 연구라 함은 자궁내 미생물총의 차이가 반복적인 유산과 조산등에 대한 원인이라는 연구가 있어 자궁내 미생물총이 host인 여성의 건강, 그리고 직접적으로는 태아의 건강과도 연관되지 않을까 하는 연구입니다.

본 리뷰지에는 10개의 자궁내 미생물총 연구 결과 논문을 소개하고 있는데 개별 연구에서 일관성있는 결과들이 나오지 않아서 본 리뷰 논문 제목처럼 자궁내 미생물총은 우리와 함께하는 파트너(거주자)인지, 어떠한 이유로 자궁으로 들어와서 잠시 머무르는 존재(여행자)인지, 악의적인 목적을 가진 존재(침입자)인지 질문을 하고 있습니다.

일관성이 보이지 않는 개별 연구의 내용들은 IVF 시술을 받은 여성들 사이에서 차이를 본 연구, 별도의 자궁이상이 없다고 판단되나 반복적인 유산 및 착상 실패를 보이는 여성들의 미생물총 차이, 자궁내막증이나 기타 다른 질환으로 자궁 수술을 받아 건강한 삶을 살고 있는 여성과 자궁내막증과 같은 질환을 가지고 있으나 아직 수술전인 여성들과의 미생물총 차이 등등등... 

얼핏보면 좀 차이가 있어보일만한 그룹을 비교를 했는데도 왜 연구마다 결과가 상이 혹은 반대의 결과를 보이는것은 실험방법의 한계, 즉 피험자의 코호트 크기(피험자 수, 논문을 보시면 아시겠지만 30명전후라서..), 시퀀싱의 한계, 먹거리와 같은 차이로 인한 환경의 차이, 샘플링방법... (물론 멋지다고 생각한 연구 디자인도 빠질수 없겠죠)

여튼 그래서 논문저자들은 각 연구의 결과들을 취합하여 다음과 같은 자궁내 미생물에 대해서 설명하고 있습니다.

1. 파트너(거주자)로써의 자궁내 미생물의 역할

 호르몬에따라(월경주기/폐경기전후) 비슷한 미생물 프로파일을 보이는 경향이 존재.

2. 잠시 거쳐가거나(여행자) 기회를 봐서 문제를 일으키는 존재(틈을봐서 돌변하는 침입자)

 수정을 위해 잠시 들어오게되는 미생물들이 있을 수 있음을 확인. 그러나 일부 인과관계는 확인되지 않았으나. 자궁 내막증을 가진 여성의 자궁 미생물 프로파일은 그렇지 않은 여성의 자궁 미생물 프로파일과 차이를 보였음.

3. 이도저도아니고 걍 백해무익한 존재(침입자)

 특정 질병을 가지고 있거나 특정 문제가 있을때 자궁/질에서 우점하고 있는것처럼 보이는 미생물 프로파일 결과가 있음.

결과는 그렇다고 쳐도 다양한 방법과 피험자들을 대상으로 pilot연구가 되어 있으니

우리는 이것을 보고 한걸음 더 나아가면 되지 않을까 합니다.

그리고 gut-brain을 시작으로 다양한 axis가 있는데 이번기회에 gut-organome 중 하나로 gut-vagina axis을 주창하는 기회를...

(와.. 역시 나와있네요.... Female Gut and Genital Tract Microbiota-Induced Crosstalk and Differential Effects of Short-Chain Fatty Acids on Immune Sequelae )

부록으로 논문에서 언급된 논문들을 SRA에서 확인해 봤는데, 혹시 사용하시려면 한번더 확인해보고 사용하시면 되겠습니다.

논문	SRA번호
Fang et al. (49)	PRJEB9626
Moreno et al. (47)	PRJNA329174
Walther-Antonio et al. (35)	SRP064295
Chen et al. (27)	PRJEB16013, PRJEB21098

저자 뒤에 표시된 괄호안 숫자는 논문내 에서 표시되는 reference 번호입니다.

출처: @sana_twice.09

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

현재 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

qiime2 -시작하기1-

qiime을 한번 실행시켜 봤습니다.

저는 완벽한 ubuntu 시스템이 아닌 윈도우에서 지원하는 linux subsystem에서 ubuntu 18.04 LTS를 설치해서 사용하고 있습니다.

그래서 약간 이글을 보시는 분들과 환경이 다소 차이가 날 수 있는걸 감안하시고 봐주시기 바랍니다. :)
-근데 지금까지 크게 차이 안나는거 같습니다. ㅎㅎ 소스 컴파일을 빡시게 해봐야 좀 차이가 나는걸 느낄수도...

WSL의 ubuntu를 시작하면 걍 shell창이 뜹니다.

그러면 qiime2를 원할하게 사용하기 위해서
conda env를 활성화 시켜보겠습니다. (물론 꼭 conda를 활성화해서 안해도 사용 가능핧!!..)

$conda activate qiime2-2019.10
하면 다음과 같이 conda env환경이 활성화 됩니다.
(qiime2-2019.10)$

다음과 같이 qiime 명령어를 실행시켜주면
(qiime2-2019.10)$qiime






















사용가능한 명령어들을 표시됩니다.
사용 방법은 모 알고 있다 싶이 qiime뒤에 명령어 입력해서 엔터치시면
(qiime2-2019.10)$qiime [command]

또 각 명령어에서 사용 가능한 sub command들을 볼 수 있습니다.

이제 qiime2를 본격적으로 사용하시기 전에 qiime glossary를 잘 일독하시고 사용하시면되겠습니다.

튜토리얼 페이지 상단에 이렇게 적혀있습니다.
Abandon all hope, ye who enter. If you have not read the glossary.
(신곡을... ㅋㅋ 센스가..) 꼭 읽으십시요!!

여튼.. 저는 glossary없이 qiime2를 접했을 때 첫 관문이 artifact, artifact를 이런 의미로  사용할 줄이야.. artifact는 원래 제거 대상아니었나요 Orz

qiime1으로 파이프라인 셋업해본 분이시라면 이건 모 신세계 아니겠습니까!!
일목요연한 명령어 리스트를 명령어 하나로 볼 수 있다니.. 아놔

그럼 다음에 또 명령어 몇개 실행해보고 다시 와보겠습니다. :)

출처: @sana_twice.09

※ 짧게라도 정리해서 일주일에 한개는 포스팅 할 줄 알았는데 이거 모 보름에 하나 겨우겨우 하네요.. ㅋ 정신을 가다듬을 정신이 없어서... 일단 되는대로 한번 업데이트 하는걸로...

qiime2 설치편

지난번 qiime2 파헤치기에 이어
이번에는 설치편!!
WSL을 이용하여 우분투 18.04 LTS에 설치하는것처럼 흉내를 내보겠습니다.

윈도우 10에서 WSL을 활성화시키고 Ubuntu를 설치하면 되겠습니다.
자세한 내용은 >여기< 를 참고해주시면 될것 같습니다.
(저도 이분의 블로그를 보고 WSL 셋팅하고 ubuntu 설치했습니다. :) )

위에서 우분투 설치후 apt-get update와 build-essential같은 패키지 설치해 주시면 다음 작업에 크게 문제가 생길일은 없으니 설치하라고 하시는거 다 설치해주시면됩니다만 생정보도구중에 특이하게 요구하는 패키지를 사용하는 프로그램 쓰시는 분이라면 알아서 하실만하리라고 봅니다. :)

여튼 ubuntu 업데이트까지 되고 나면

qiime2를 설치하기 위해 conda를 설치하시면됩니다.

conda는 miniconda2를 사용하였고 자세한 설치 화면은
>여기< 를 참고하시면되겠습니다.

$wget https://repo.anaconda.com/miniconda/Miniconda2-latest-Linux-x86_64.sh
$sh Miniconda2-latest-Linux-x86_64.sh

miniconda2를 설치한 이후에는 qiime2의 설치 안내 페이지에 나와 있는대로
진행해주시면됩니다. 자세한 내용은 >여기< 에 나와있습니다.

$conda update conda
$conda install wget
$wget https://data.qiime2.org/distro/core/qiime2-2019.10-py36-linux-conda.yml
$conda env create -n qiime2-2019.10 --file qiime2-2019.10-py36-linux-conda.yml

conda에 설치된 qiime2를 사용하려면 qiime2가 설치된 conda환경인 qiime2-2019.10을 활성화 시켜주면됩니다.

$conda activate qiime2-2019.10

활성화된 환경을 비활성화 하려면 다음과 같이..

$conda deactivate

일단 설치는 했습니다.

다음에 사용해보아요~ :) 제발~

출처: @sana_twice.09

Windows Subsystem Linux 입문

?? 윈도우용 리눅스 하위 시스템 설명서가 아니라
Linux용 Windows 하위 시스템 설명서라고?

번역이 다른건지 문화가 다른건지 그건 중요한게 아니고
윈도우에서 cygwin도 아니고 버추얼박스도 아닌데 리눅스 명령어를 그것도 우분투를 사용할 수 있다고?
이런 말도 안되는 일이 이미 꽤 예전부터 사용되고 있었는데
저는 이번달에 되서야 알았고..

별도 서버를 접속해야 테스트를 할 수 있었던걸..
이젠 손쉽게 할 수 있게 되었습니다.

윈도우10에서 ubuntu에 접속했습니다! (아... 전 이 화면을 좋아합니다.)

역시 Dell 사용자로 정품 윈도우를 사용하니 좋은 날이 있군요
모 어둠의 경로로 해서 사용해도 가능은 할거 같은데..

WSL도 WSL1과 2로 나눠지는거 같은데...
일단 우분투 설치도 간편하게 잘 했고
이것저것 설치 또한 간편하게 다운받아서 했고
전혀 시스템에 별도로 오버헤드 걸리는것 없이 자연스럽게 작업을 하고 있습니다.

다만 일정 작업을 하면 자원을 완벽히 공유하고 있어서
순간 빡센 작업을 시키면 서류작업같은거 전혀 안되는게 단점입니다. ㅎㅎ


그래서 조만간 지난번에 진행하지 못했던 qiime2 설치편을 이어서....

출처: @sana_twice.09

엄마한테서 니가 왜 나와

중국에서 비 침습성 산전 검사(NIPT)를 활용해서 흥미진진한 연구를 진행해서 2018년도에 투고한 논문을 한번 훑어봤습니다.

Genomic Analyses from Non-invasive Prenatal Testing Reveal Genetic Associations, Patterns of Viral Infections, and Chinese Population History
>논문은 여기로<


NGS시대에 유전학연구를 할 때 가장 문제되는게 표본 크기가 제일 문제되고 있죠, 시퀀싱을 하자니 기본적인 depth가 필요한데 적절한 depth로 하자니 비용때문에 표본 크기를 줄여야 되고.. 표본 크기를 줄이자니 그정도 시퀀싱해서 얻어진 결과가 잘 될까?

그래서 중국 연구진(물론 미국과 덴마크 연구진도 합세했습니다.)이 NIPT 데이터를 이용해서 유전체 연구를 진행해 보았답니다.

NIPT는 shallow sequencing의 대명사죠, 그럼에도 불구하고 전세계적으로 스크리닝방법으로 정착을해서 현재도 굳이 침습적 산전진단을 받지 않아도 되는 산모들에게 대안이 되어 주고 있습니다.
(오늘은 NIPT 성능과 위음성과 같은 것은 논외로 하겠습니다.)

산모들의 NIPT 데이터를 사용해서 유전학 연구를?
이게 될까? 저도 될줄은 몰랐는데 이분들 일단 했으니 논문이 나왔겠죠..

대한민국에도 엄청 많은 NIPT 데이터가 쌓여 있을텐데.. 한번 해봐도 나쁘지않을 것 같습니다. (물론 동의서가... ㄷㄷㄷ )

결론적으로 슬쩍 읽는 NIPT 가지고 거기다가 두명혹은 그 이상의 정보가 섞여 있는데 그거가지고 모 뭐할꺼야? 라고 하신다면 큰 오산!
슬쩍 읽는가지고도 집단유전학, Association mapping study 그리고 (아놔 이거 모야?인데 진짜 나온다는.. ㅋㅋ ) 인간 혈장내 바이러스 연구를 할 수 있다고 합니다.

일단 NIPT로 얻어진 allele 정보가.. 진짜 쓸만한가 알아보았슴다.
NIPT 데이터(CMDB)에서 call된 SNP들이 기존의 1KG CHN,EAS 및 전체에서 얼마나 확인되었는지, CMDB에서 novel variant들이 있는지 확인했고
non-reference allele freq를 가지고 1K genome의 CHN allele과 CMDB의 allele를 비교해보았습니다.

물론 Novel은 적게, gnomAD, 1KG 데이터와 많이 공유하고, 상관관계는 높습니다.


변이에 대한 결과를 확인 했으니 이제 중국의 유전학 썰을 한번 풀어보는걸로..
분석에 걸림돌이 되는 오류같은것들은 적절히 제외해주고 한족에 대해서 한번 확인해봤습니다.

31개 행정구역의 4만5천명 한족이라고 확인된 샘플들을 PCA분석을 해보 았는데 위도에 따라 차이가 많이 보였습니다(이 결과는 이전 연구 결과와 동일한.., 동서에 의한 차이는 별로.. ,논문내에서는 이주정책에 의해 설명되는것 같다고 하는데 정확히 무슨 말인지는.. 1996년 논문을 봐야할듯 하네요. ). 그러나 그 외 소수민족에서는 걍 경도, 위도 상관없이 다 차이가 있었다고 하네요. 그래야 소수민족이지 않을까하는 생각이기도 합니다.
가장 차이가 나는 집단은 신장지방의 Kazakh, Uyghur, 내몽골 지역의 Mongol이었습니다.
(일단 언어가 다르다고 하네요)

그리고 키와 BMI와 같은 Multiple Complex Traits의 Phenotype - Genotype Association을 확인해봤슴다.
지금까지 연구된것들과 비슷한 SNP들이 확인되었고 CBK, GIANT, UK BioBank 대규모 set들과 비교해봤을 때도 크게 다른 값을 보이지 않았습니다. 그리고 NIPT라는 특별한 샘플이라는 점을 고안하여 쌍태아와 관련된것도 한번 확인해 보았답니다.
(산모 연령에 대한 것도 확인했다고 했는데... 이건 유전적이라기 보다는... 사회학적인 부분이지 않을까했는데.. 모 연관성이 높은 SNP를 찾긴 찾았는데 저는 패스 하도록 하겠습니다.)
중국의 NIPT 샘플중 470여명의 쌍둥이들이 있었는데(초음파로 확인) NRG1의 rs12056727d이 꽤나 흥미로운 좌위로 확인되었습니다. 이 SNP는 UK BioBank에서 갑상선 항진증과 관련이 있는것으로 확인되었습니다. 그리고 갑산성 기능은 불임과 관련이 있습니다. NRG1유전자에 대한 내용은 일단 넘어가더라도 쌍둥이를 임신한 산모는 갑상선 자극 호르몬 수준이 낮아지는 경향이 있습니다.


그리고 오늘 포스팅 제목에 부합하는 내용인, 산모 혈장에서 보이는 circulating viral DNA!
왜왜왜 나오냐고

일단 NIPT에서 수행된 시퀀싱 결과에서 human에 mapping되지 않는 read들을 따로 모아서 NCBI 바이러스 DB에 확인해봤고 그 결과를 Figure5의 A,B에 정리했습니다. 그렇지만 우리가 잘못 분류하거나 샘플이 오염 되었을 수도 있기에 검토하였다. 대부분은 특정 지역에 국한되었고, 이미 인체내 genome에 내재화 되어 있는 HCV서열이 있는것을 확인했지만... 그 대상 숫자가.. 고작 3명..

근데 중국에서만 그런건 아니고 유럽쪽 연구에서도 바이러스가 나왔었다능
-물론 바이러스 종류가 좀 다른 것으로보입니다.
실험 디자인이 달라서 그런것인지 진짜 지리적으로 차이가 있는것인지 모 좀더 확인해봐야 하겠지만... 크.. 인구도 많고 땅덩어리도 작지 않은 중국 아니겠습니까.
HBV관련해서 민감도도 계산했고 지리적 분포를 확인해서 지리적으로 확인되는 바이러스 종류들이 다르다는 것을 확인했습니다.

그래서 NIPT 데이터를 잘 활용한다면 공중보건에도 꽤나 유용할것 같다는게
저자들의 논지였고 현재 NIPT 서비스를 하고 있는(그러나... 내년에도 하고 있을거라는 보장은... 어느 이유로든지 밝지 않은) 본인에게도 걍 헛소리로 보이지는 않고 진짜 유용할것같다능...

ps. 근데 일주일에 논문하나 읽기가 이렇게 어려웠나;;;

출처: @sana_twice.09

저쪽집이 좋지만 우리집도 잘해요, 시퀀싱

가능하면 일주일에 하나씩은 업데이트 하려고 했는데 여윽시..
그건 어려운것 같네요 ㅎㅎ
그래도 되는대로 논문읽고 일주일에 한번씩 업데이트 하는걸로 :)


오늘은 금년 6월달에 남중국과학대학에서 Scientific reports에 투고한 논문 되겠습니다.

제목은 Systematic comparison of germline variant calling pipelines cross multiple next-generation sequencers >여기를 방문하세요<

현존하는 1빠 시퀀서대비 가격도 저렴 시약도 저렴한  BGISEQ과 MGISEQ 성능 비교 테스트인데 결론은 왜 Strelka2가 적절한 분석 파이프라인으로 권장한다인지..

-일단 BGISEQ과 MGISEQ의 라이브러리 제작 및 시퀀서 방식을 뒤로하고 성능만 봅니다.
덤으로 Tianhe-2라는 슈퍼컴퓨터 자랑도 -


3개 콜러(GATK4, Strelka2, Samtools-Varscan2)를 가지고 WES, WGS를 비교해보니 WES데이터는 시퀀서 및 콜러별로 높은 일치성을 보이는 반면 WGS는 그러지 못했습니다(니들도 WGS가 WES처럼 높은 일치성을 보일거라고 생각안했잖아ㅋㅋ 어디서 약을.. 그래도 논문은 나왔기에 괜찮습니다 Orz )

Figure 1. 우리는 여러분이 가장 많이 사용하고 권장하는 콜러를 지금 있는 그대로 분석을 돌려봤습니다.

Sequencing Samples	Bases(Gbp)	Read(x106)	Clean rare	>Q20	>Q30	GC	Mean coverage
BGISEQ500-WES	29.41	294.3	0.41%	96.72%	89.14%	49.75%	328.49X
MGISEQ2000-WES	16.34	163.55	0.25%	98.18%	92.08%	49.71%	129.40X
HiSeq4000-WES	41.93	283.7	4.46%	97.36%	93.01%	50.63%	395.17X
NovaSeq-WES	25.88	178.87	2.25%	95.33%	92.67%	49.73%	241.52X
BGISEQ500-WGS	126.86	1270.02	1.76%	93.73%	83.33%	41.76%	41.03X
MGISEQ2000-WGS	137.36	1374.87	0.21%	96.17%	88.19%	41.76%	45.13X
HiSeq4000-WGS	191	1276.1	8.25%	95.90%	90.11%	41.69%	58.00X
NovaSeq-WGS	98.3	657.45	1.28%	95.89%	93.86%	41.61%	28.96X
HiSeq Xten-WGS	134	894.58	7.29%	94.50%	87.63%	40.71%	38.93X

Table 1. 우리 필요한 만큼 시퀀싱 잘 했어요

Figure 4,5 WES관련 작업 시간 및 결과 정리
Figure 6,7 WGS관련 작업 시간 및 결과 정리
(이미지 생략)

그래서 WES와 WGS를 각각 콜러의 조합에서 분석한 결과 SNP는 일관적으로 잘 call하였고 InDel은 일관적이지 못했다.
플랫폼별로 보면 SNP는 BGI플랫폼이 InDel은 일루미나 플랫폼이 더 나았다.
이거슨 시퀀싱할때 read 길이를 BGI 플랫폼과 illumina 플랫폼의 길이를 각각 100PE, 150PE로 해서 그렇다는 이유를... (그럼 왜 BGI플랫폼은 150PE로 안하고..??)

그리고 시퀀싱 뎁스 운운하는데.. 결론은 추가테스트 및 다른 분들이 더 해줬으면 하는걸로..
그리고 각 플랫폼에서 분석 툴의 성능 비교는 Strelka2가 다른 2개 분석 방법도나 나은걸로
순위를 따지자면 Strelka2 > GATK > Samtools-VarScan(SV) (모 다들 예상하셨다 싶이..)
InDel을 call하는 결과가 좀 차이가 있었는데 NovaSeq-SV에서 23개의 small variant를 call했는데반해 X-Ten-SV에서는861개의 small variant을 call했.. (BGI플랫폼은 갑자기 사라지고..)
그리고 마무리는 germline의 SNP, InDel call 능력은 높은 일치성이 있는 걸로 마무으리~

종합적으로 Strelka2가 최적의 분석 파이프 라인
응? 이거 시퀀서 비교 아니었어?
응 아니야, Strelka2 좋아요 꾹! 구독아니 github 꾹!
결국 이렇다고합니다.

이 논문의 의의는 테스트 해볼 비교 set이 생겼다는것에 ...  :)

출처: @sana_twice.09

산모의 건강상태가 NIPS에 미치는 영향

Noninvasive prenatal testing for aneuploidy using cell-free DNA – New implications for maternal health

간만에 하는 업무와 갱장히 밀접한 글 하나 투척이라기 보다는
간단히 정리 정도...


NIPS/NIPT, 비침습적 산전 스크리닝 혹은 테스트의 약자로 산모의 혈액내에 있는 태아의 cfDNA를 사용하여 태아에 문제가 있는지 확인 할 수 있는 기술로써 홍콩대학의 Denis Lo 옹께서 산모 혈액내에 태아것도 있다. 그리고 산모와 태아의 것을 구분 할 수 있다고 밝혀주셔서 전세계적으로 빠르게 확산되었죠.


지속적인 연구를 통해 Trisomy 13/18/21에 대해서 우수한 정확도를 보여주고 있습니다만, NIPS/NIPT를 하다보면 이게 검사가 잘 안되는 경우가 있습니다.

그러나 어떠한 이유인지 분석 할 수 없는 경우나 알 수없가 없었습니다.
그래서 지속적으로 축적된 임상 경험을 바탕으로 NIPS/NIPT와 산모의 건강과 관련하여 상호관계가 있을것으로 추측하였고 그래서 한번 어떤 것들이 NIPS/NIPT의 성능에 영향을 미치는지 조사해본 논문하나를 살펴보았습니다. :)


산모의 비만 (Obesity)

NIPS/NIPT를 할 때 의뢰서에 키와 몸무게 란이 있습니다. 이는 BMI이 치수로 산모의 체질량을 확인하는데 BMI가 높을수록 NIPS/NIPT에 필요한 태아의 cfDNA의 비율이 떨어지기 떄문에 NIPS/NIPT가 실패할 확률이 높아진다고 보고되고 있습니다.
지방세포에서 세포 사멸과 괴사로 혈중 내 태아의 cfDNA보다 산모의 cfDNA 양이 상승하여 태아의 cfDNA 비율을 떨어트린다는 것입니다.


자가면역질환 (Autoimmune disease)

자가면역질환은 문자 그대로 자기 자신을 공격하는 질환으로 세포 주기가 굉장히 짧아지는것이 문제가 되어 임신하지 않은 전신성 홍반성 루프스(SLE) 환자에서 혈중 cfDNA가 증가하는 것을 확인했습니다. 그래서 SLE가 있는 임산부의 NIPS/NIPT결과를 해석할 때는 주의를 기울여야 한다고 말하기도 하였습니다.


태아 분획에 영향을 줄 수 있는 것들 (Other maternal influences on fetal fraction)

표2(논문 참조)에 나와있듯이 태아 분획과 다양한 산모의 특징(인종, 흡연여부, 고혈압, 나이등등등)들은 통계적으로 유의한 연관성을 가지고 있으나 이러한 산모의 특징들은 상반된 결과를 나타내는 결과들도 있기 때문에 덜 중요해 보인다고 합니다.
(인종 ㅂㄷㅂㄷ, ㅂㄷㅂㄷ 했다고 인종차별자는 아니니 오해 없으시기 바랍니다.)


산모가 받은 처방 (Maternal medications)

산모가 복용하는 약 또한 NIPS/NIPT의 결과를 방해 할 수 있습니다.
LMWH(Low molecular weight heparin)은 NIPS/NIPT에 부정적인 영향을 미치는 첫 물질이었습니다. heparin과 NIPS/NIPT의 부정적인 결과에 대해서 정확한 메카니즘은 알려져 있지 않지만 NIPS/NIPT 검사를 위해서는 heparin 복용을 일정 시간 중단해야 합니다. 그리고 IVIG(Intravenous immunoglobulin)의 경우 cfDNA가 증가하는 것처럼 보고 되었지만 NIPT의 분석 성능에는 확인되지 않았습니다.


장기 이식 받은 산모 (Organ transplant recipients)

장기 이식을 받은 산모의들의 NIPS/NIPT는 고려를 해봐야 합니다. 특히나 장기 이식받은 장기가 남성 기증자라면 특히나 주의해야한다. 남성 기증자의 장기는 Y염색체 특이적 cfDNA서열을 만들어내어 잘못된 판단을 유발 할 수도 있다. 그리고 NIPS/NIPT는 염색체 이상이 있는 산모가 하면 정확하지않은 결과를 받아볼 수 있기 때문에 NIPS/NIPT 전 상담을 꼭 진행해야 합니다.


※ 정정 :하단은 NIPS/NIPT를 해서 산모로부터 우연히 알게되는 것에대한 내용입니다.

산모의 성염색체 상태 (Maternal sex chromosome conditions)

NIPS/NIPT에서 X염색체의 정확도는 13/18/21과 같은 상염색체보다 본질적으로 좋지 않습니다. 다른 이수성에 비해 X 염색체 이상에 대해 더 높은 태아 및 태반 모자이크 현상이 있습니다. 산모에 따라 나이와 연관있는 X 염색체 손실이 일어나기도 하고 일부 산모에게서는 정상 XX에서 저수준으로 XO/XX 모자이크로 바뀌기도 합니다.또한 많은 SCA(Sex chromosome aneuploidies)는 불확실한 NIPS/NIPT결과를 초래하기도 합니다. NIPS/NIPT에 의해 의심된 181개 SCA case중 16 case (8.5%)가 산모에 의해 기인한 것으로 확인되었다고 합니다.


산모의 상염색체 이상 (Maternal autosomal abnormalities)

산모의 상염색체이상도 잘못된 NIPS/NIPT 결과를 발생시킵니다. 이런 경우 산모가 저수준의 모자이크 T18이나 염색체 13 혹은 18번에 부분적으로 중복이 있어 발생할 수 있습니다. 최근 NIPS/NIPT 분석 중에 22q11.2와 같은 미세결실 증후군을 탐지하고 있는데 한 연구에서 고위험의 미세 결실 결과 중 60%(20 case/32 case)가 산모에 의한 것이지 않을까 의심되고 있습니다. 그러나 다른 NIPS/NIPT 제공 업체의 경우 고위험 미세 결실 결과중 2%만이 산모에 의한 것으로 확인(?이 진짜 된건지 산모를 검사를 했는지는 까지는)됐다고 합니다. 그런데 이런 차이가 나는 이유는 NIPS/NIPT분석 방법이나 reference의 편향성이나 인구 특징들에 따라 달라지기 때문에 분석 결과에 대한 비교까지는 어려운것 같다고 합니다.


산모의 종양 (Maternal malignancy)

종양세포는 NIPS/NIPT 분석에 사용하는 혈장내의 cfDNA를 기여하고 있다는 것에 누구도 부정은 못할 것 입니다. 종양이 NIPS/NIPT 결과에 영향을 준다는 사례 보고는 여기에서 설명되었습니다. 125,426명의 산모 집단에서 산모의 종양으로 인해 NIPS/NIPT결과가 불일치한 산모 10명에 대해서 정보를 공개했습니다. 산모들의 종양/암 유형은 lymphoma, leukaemia, colorectal and anal cancers 이었다고 합니다.
그리고 산모들은 추가적인 생물정보학적 분석으로 비특이적인 CNV들을 보여줬는데 NIPS/NIPT에서 관찰된 비특이적인 패턴들은 종양/암 치료후 사라져서 종양/암과 NIPS/NIPT와의 관계를 확인 할 수 있었습니다.


그래서 종합적으로 요약하자면 임상적으로 NIPS/NIPT에 대한 경험과 자료가 증가함에 따라 산모의 건강과 밀접한 관련이 있다는 것을 확인 할 수 있었고 비만, 면역질환 및 LMWH 치료와 같은 NIPS/NIPT 결과에 영향을 줄 수 있는 산모의 상태에 대해서 인지하고 있어야합니다(안그러면 위양성 결과를 받아들고 옳지못한...).
NIPS/NIPT결과로 얻어진 2차 산모의 위험 소견에 대해 산모에게 알리는 것도 중요한 윤리적 고려사항이라고 생각됩니다.


이렇다고합니다.

간만에 업무와 관련되 내용 포스트 좀 해봤습니다.

출처: sana_twice.09

bam파일을 가장 쌈빡하게 잘 다루는 방법

NGS 시대를 살면서 bam파일은 어떻게 땔래야 땔수없는
껌딱지같은 녀석이 되어버렸습니다.

다양한 정보를 가지고있기는 한데 그 안에서 정보를 추출같은 핸들링하기가 쉽지 않은게 사실이긴 합니다.

그래서 여기 sam파일 bam파일을 samtools보다 빠르게  처리해주는 툴하나
가지고 와봤습니다.

모 이미 samtools대신에 잘 사용하고 계실겁니다.

sambamba라고 (현재 버전은 v0.7.0 입니다. 논문은 여기에 있습니다.)

sambamba는 samtools와 비슷한 대부분의 기능을 가지고 있습니다.
view, sort, index, merge, slice, flagstat, markdup
그리고 사용법도 모 크게 어렵지 않습니다. :)

sambamba documentation을 보시면 모 크게 어렵지 않게 손쉽게 사용하실 수 있을겁니다. :)

출처: sana_twice.09

Trimmomatic 설치 및 사용

Trimmomatic Site

Trimmomatic Manual

Trimmomatic은 cutadapt와 함께 Illumina 시퀀서들의 adapter들을 제거하고 trimming하는데 널리사용되는 tool되겠습니다.

설치방법은..
>wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.36.zip
>unzip Trimmomatic-0.36.zip
Trimmomatic-0.36폴더 밑에 파일 생성
>cd Trimmomatic-0.36
>java -jar trimmomatic-0.36.jar -h

참 쉽죠!!
물론 java가 설치되어 있어야한다는게 함정


Quick Start

Paired End일때
java -jar trimmomatic-0.35.jar PE -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 ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36


Single End일때
java -jar trimmomatic-0.35.jar SE -phred33 input.fq.gz output.fq.gz ILLUMINACLIP:TruSeq3-SE:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

위의 Quick Start에서 사용한 내용은 paired와 single 차이일뿐 옵션은 동일합니다.
찬찬히 뜯어보면 다음과 같다고 하네요

ILLUMINACLIP:TruSeq3-PE.fa:2:30:10
- 일루미나 아답타 서열을 제거하는 옵션입니다. 자세한 내용은 하단에..
LEADING:3
- 서열의 앞쪽을 기준으로 quality (여기서는 3) 이하의 N개 서열을 제거합니다.
TRAILING:3
- 서열의 뒤쪽을 기준으로 설정한 quality (여기도 3) 이하의 N개 서열을 제거합니다.
SLIDINGWINDOW:4:15
- 4base씩 확인하면서 평균 quality가 설정된 기준 quality (여기서는 15)보다 이하일 경우 제거합니다.
MINLEN:36
- 위의 단계들을 거치면서도 살아 남은 서열 길이가 기준 길이 (여기서는 36bp)보다 짧으면 (이하) 아예 read를 삭제합니다.



사용 가능한 옵션

ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>
- fastaWithAdaptersEtc: trimmomatics에서 TruSeq2와 TruSeq3서열은 제공하고 있고 그외에 아답터 서열을 사용하고 싶으면 fa파일을 넣어주면 됨
- seed mismatches: 허용가능한seed 서열의 mismatch 개수
- palindrome clip threshold: palindrome trimming방법은 Paired End로 시퀀싱을 했을 때 적용할 수 있는 방법으로 아답타 서열이 traget read 앞뒤에 존재해서 첫번째 read와 두번째 read에서 역 상보적으로 동일한 read들이 발견되는 경우 해당 서열들을 삭제하는 방법입니다. (참고 palindrome란? )
- simple clip threshold: 설정된 기준 값에 부합하는 충분히 정확한 match가 확인되면 적절히 clipping함

SLIDINGWINDOW:<windowSize>:<requiredQuality>
설정된 windowSize 서열들의 평균 qaulity가 requiredQuality 값보다 작아지면 작아진 서열들 이후를 제거합니다.

LEADING:<quality>
read 앞쪽을 기준으로 설정된 quality보다 낮은 quality를 가진 서열들을 제거합니다.

TRAILING:<quality>
read 뒤쪽을 기준으로 설정된 quality보다 낮은 quality를 가진 서열들을 제거합니다.

CROP:<length>
read 앞쪽을 기준으로 설정된 length만큼 서열들을 보존한다.

HEADCROP:<length>
read 앞쪽을 기준으로 설정된 length만큼 서열을 삭제한다.

MINLEN:<length>

TOPHRED33 quality를 phread 33으로 변환

TOPHRED64 quality를 phread 64으로 변환


그리고 중요한 점 하나! trimmomatic는 옵션에 순서가 있다는 사실! 저도 아직 테스트 해보지는 못했는데 옵션이 순서대로 적용 된다고 합니다. 설명서에 그렇게 써 있으니 맞겠죠. 그래서 일단 어답터 서열을 먼저 제거하라고 합니다. :)

그럼 즐거운 trimming 작업 되시길..

@sana_twice.09