Be great

tqdm을 사용해서 얼마나 진행되었는지 작업 진행 정도를 표시함.

전체 양을 알 때와 모를 때를 나눠서 표시 할 수 있음.

전체 양을 모를 때 

from tqdm import tqdm
with open(filename) as f :
       for index, line in enumerate(tqdm(f, unit='reads', unit_scale=True, mininterval=1)):
               continue

결과물 :

58.9Mreads [00:24, 2.38Mreads/s]

3회 평균 소요 시간 21.7초

전체 양을 알 때

from tqdm import tqdm
with open(filename) as f :
        lines = f.readlines()
        for index, line in enumerate(tqdm(lines, total=len(lines), unit='reads', unit_scale=True, mininterval=1)):
                continue

결과물 :

100%|██████████████████████████████████████████████████████████| 58.9M/58.9M [00:15<00:00, 3.88Mreads/s]

3회 평균 소요 시간 30.6초 

f.readlines() 함수로 파일 전체를 읽어서 사이즈를 계산하고 진행 했을 때는 이미 메모리에 내용이 올라왔기 때문에 시간 당 읽는 줄 수는 빠르지만 파일을 읽는데 드는 시간으로 인해 총 시간은 더 느림 하지만 전체 진행율을 알 수 있다는 장점이 있음.

conda 설치 후 채널을 추가해야 최신 버전의 프로그램을 갱신해서 쓸 수 있다.

conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge

import pandas as pd

df = pd.DataFrame()

#make dataframe from dictionary
tmp_df = pd.DataFrame([foo_dic], index=id)

#sum data by raw
total_count = tmp_df.sum(axis=1)[0]
tmp_df = tmp_df.div(total_count)

#concat multiple dataframe parellel
concat_df = pd.concat([df1,df2],axis=1).fillna(0)

#specific columns contain letter 'test'
df = df.loc[:,df.columns.str.contains('test', regex=True)]

#merge dataframe consider index
merge_df = pd.merge(df1, df2, left_index=True, right_index=True)

log 파일 작성 모듈 logging.

import logging

logger = logging.getLogger()
logger.setLevel(logging.DEBUG)
formatter = logging.Formatter('[%(levelname)s %(asctime)s] %(message)s',"%Y-%m-%d %H:%M:%S")

stream_handler = logging.StreamHandler()
stream_handler.setFormatter(formatter)
stream_handler.setLevel(logging.INFO)
logger.addHandler(stream_handler)

file_handler = logging.FileHandler(f'my.log')
file_handler.setFormatter(formatter)
file_handler.setLevel(logging.DEBUG)
logger.addHandler(file_handler)

logging.info(f'Read Database File')
logging.debug(f'Read Database File')

handler를 여러 개 만들어서 하나는 stdout 다른 하나는 my.log 파일로 만들고 level에 따라 출력 범위를 다르게 조절한다.

위의 예시에서는 'DEBUG', 'INFO', 'WARNING', 'ERROR', 'CRITICAL', 중에 debug는 파일로만 생성되도록 설정되었다.

개발 단계에서는 stream_handler를 DEBUG로 놓고 진행하다가 개발 완료시 INFO로 수정하면 원하는 부분만 출력하도록 조정 가능하다.

a = subprocess.Popen(cmd, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT).communicate()[0].decode('UTF-8')
logging.debug(f'{cmd}\n{a}')

subprocess와 연결해서 command를 입력하고 나오는 출력물까지 debug로 한 번에 연결 할 수 있다.

read 서열에서 error를 1이하로 허용하는 내에 BESTMATCH를 찾아 시작과 종료지점 그리고 매치되는 서열을 확인하는 코드. error는 mismatch, insertion, deletion을 의미한다.

import regex

primer_seq, read_seq

regex_primer_seq = fr'({primer_seq}{{e<=1}})'
match_object = regex.search(regex_primer_seq, read_seq, regex.BESTMATCH)

match_start, match_end = match_object.span()
match_seq = match_object.captures()

anaconda의 버전 업데이트 후 아래와 같은 에러 발생

Collecting package metadata (current_repodata.json): failed

ProxyError: Conda cannot proceed due to an error in your proxy configuration.
Check for typos and other configuration errors in any '.netrc' file in your home directory,
any environment variables ending in '_PROXY', and any other system-wide proxy
configuration settings.

https://github.com/conda/conda/issues/9497 페이지에서 확인해보니 urllib3 의 버전 변경으로 인한 문제로 보임.

conda config --set ssl_verify False

위 명령어를 작성하면 .condarc 파일이 home directory 밑에 생성되는데 파일을 아래와 같이 수정한다

ssl_verify: False
proxy_servers:
  https: http://[PROXY]:[PORT]/

https: https: 에서 https: http로의 변경.

NGS 데이터에서 adapter 서열과 primer 서열에 따른 read 분류.

아래는 pseudo code 이므로 적절한 변환이 필요하다.

import regex

#ambiguous base list
ambiguous_base_dic = {"N":"ATGC","R":"AG","Y":"TC","K":"GT","M":"AC","S":"GC","W":"AT","B":"CGT","D":"AGT","H":"ACT","V":"ACG"}

#if ambiguous base occured, change to regular expression format
for word, initial in ambiguous_dic.items():
	primer_f = primer_f.replace(word, "["+initial+"]")
    
#index sequences are allowed one mismatch
primer = r"(^{0})".format(primer_f)+"{e<=1}"

#if primer sequence is inside sequence, print out.
if regex.findall(primer, sequence):
	print(sequence)

Kaggle이란 데이터 과학자와 머신 러닝 사용자의 스킬과 정보 공유를 위해 만들어진 커뮤니티이다.

주기적으로 열리는 Competetion을 통해 상금과 명성을 얻을 수도 있고 

Courses 에서 학습을 할 수도 있다.

상당히 공신력 있는 커뮤니티라서 입상 실적은 ML/AL 경력으로 인정해주는 회사도 있다.

Kaggle korea 페이스북 커뮤니티도 있으니 참고하면 좋을 듯 하다.

https://www.facebook.com/groups/KaggleKoreaOpenGroup/

임산부의 플라즈마에서 태아의 cell-free DNA의 존재가 밝혀지면서 NIPT (Noninvasive prenatal test)가 만들어졌고 태아의 비이상성을 검출하기 위한 방법으로 널리 사용되고 있다.

하지만 임산부의 플라즈마에서 태아의 DNA 비율이 어느정도인지에 따라 검사의 정확도가 좌우되기때문에 정확한 비율을 알아내는 것이 필요해졌고 여러 종류의 생물정보학 분석 방법이 사용되었다.

1. Y 염색체 기반

예전에는 SRY, DYS14, ZFY 같은 유전자에 마커를 사용하여 측정하였으나 최근에는 Y 염색체 전체의 비율을 계산한 것으로 FF로 계산한다. 간단하고 정확한 편이나 남아에게만 적용할 수 있는 한계가 있다.

2. 부계 Genotype 정보를 이용

SNP 중에 부계에서 A/A, 모계에 C/C 처럼 homo로 가지고 있는 경우 태아는 A/C hetero를 가질 수 밖에 없다. 전체 서열 중에 A가 가지는 비율을 계산하면 FF를 추정 할 수 있다. 이 방법도 정확도가 높은 편이지만 부계의 SNP정보를 필요로 하다는 점에서 추가 비용이 필요하다.

3. high-depth 시퀀싱 데이터

부계 genotype 정보가 없어도 충분한 depth가 확보된다면 나오는 결과를 기반으로 모계 genotype과 태아의 genotype을 추정 할 수 있다. 대략 ~120x 정도의 depth일 때 신뢰할 만하다는 연구 결과가 있다.

4. 모계 genotype 정보를 이용한 shallow-depth 시퀀싱

모계 SNP 정보를 정확히 알고 있다면 depth가 많지 않아도 FF 예측이 가능하다. 모계에서 homozygous site를 추려낸 후 이 중에 hetero 서열을 가지는 곳의 비율을 계산하면 된다. 2번과 비슷하게 SNP 정보를 microarray 등을 통해 미리 확보해야 하는 단점이 있다.

5. 시퀀싱 데이터를 기반으로한 shallow-depth 시퀀싱

50KB 단위로 염색체의 영역을 나눈 뒤 read count를 계산하여 FF값을 추정하는 SeqFF 라고 불리우는 방법이 있다. Enet과 reduced-rank regression 모델과 coefficient를 계산했을 때 각각 0.932와 0.938의 r값을 가진다. 이 방법은 실제 FF가 5% 미만일 때는 정확도를 보장할 수 없으며 대용량 데이터로 예측 모델을 만들어야하는 단점이 있다.

6. 메틸화 마커

포유류는 CpG 염기에서 C가 ~70% 정도 메틸화 되어 있다. 메틸화가 되는 위치에 따라 세포의 역할이 바뀌는데 태반 특이적 메틸화 마커인 RASSF1A 프로모터 서열 등에서 태아에서만 메틸화 되는 영역을 찾아, 모계 DNA 서열에서의 메틸화 되지 않는 영역과 비교하여 FF를 계산한다. Y 염색체 기반 FF 예측과 0.85의 correlation을 가지는 것을 확인하였으며 대규모 데이터셋에서 깊이있는 연구가 필요하다.

7. Cell-free DNA 사이즈 기반

태아의 DNA는 엄마의 DNA보다 짧은 편이다. 따라서 FF 값이 높다면 짧은 DNA의 비율이 늘어날 것이다. paired-end 데이터에서 100-150bp와 163-169bp의 DNA 사이즈 구분을 통해 Y염색체 기반과 0.827의 r값을 가지는 결과를 만들었다.

8. Cell-free DNA 뉴클레오좀 기반 

Reference - 

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5343987/

Noninvasive prenatal testing (NIPT)는 태아의 유전적 이상을 찾아내는 방법 중 하나로 임부의 피에 떠다니는 cell free DNA(cfDNA)를 검사한다. 침습적 방법의 경우 주사를 직접 태반에 찔러 넣는 방식으로 낮은 확률의 사산이 생길 수 있는것에 비해 안정적이지만 검사의 정확도가 상대적으로 낮은 단점이 있다.

cfDNA는 약 200bp의 길이로 일반적으로는 세포가 죽을 때 생기는 부산물이며 혈류를 따라 이동한다. 임신 중에는 태반 세포의 일부도 cfDNA가 되는데 태반의 DNA는 태아의 DNA와 일치하기 때문에 이 DNA를 NGS 방식으로 시퀀싱 하여 태아의 DNA를 분석할 수 있다. 

NIPT는 주로 aneuploidy 같은 염색체 이상 질환에 대해 검사한다. 특히나 주로 보는 것은 trisomy 21(다운 증후군), trisomy 18, trisomy 13 그리고 X, Y 염색체의 수 이상 등이 있다.

염색체 특정 영역의 deletion 또는 duplication도 확인 하기도 한다.

엄마의 피에서 태아의 DNA가 차지하는 비율을 "fetal fraction" 이라고 부른다. 10주 이상일 때 보통 4퍼센트 이상이며 만약 이보다 낮다면 검사의 정확도가 낮아져 위양성 결과가 나올 수 있다. 샘플링 에러, 엄마의 비만 또는 태아의 비정상적인 부분이 있다면 fetal fraction이 낮게 나올 수 있다.

염색체 이상을 찾아내는 방법은 태아와 엄마의 모든 cfDNA 조각의 개수를 세고 모든 염색체의 비율이 비슷하다면 염색체 이상이 있을 확률은 낮은 것으로 본다. 만약 특정 염색체의 조각이 많다면 trisomy를 의심해 볼 수 있다. NIPT는 정확도가 높지 않기 때문에 추가적인 검사로 양성 판단이 필요하다.

Reference -

https://medlineplus.gov/genetics/understanding/testing/nipt/