역시나 회사명은 밝히지 않음.
어제 면접보면서 코테 얘기가 나와서 블로그에 올려도 되겠느냐고 여쭤봤고, 그 쪽이 커리어에 도움이 될 테니 올려도 된다+올리면 안 되는 거면 따로 언질을 주겠다는 답변을 받았음.
문제가 이거 말고 하나 더 있는데 그거는 바이오파이썬을 쓸 일이 없고(이것도 안썼지만...), 뇌에 블루스크린이 떠서 이해하느라 시간이 좀 걸려서 아마 구현은 빨라야 내일... 잘못하면 다음주까지도 갈 것 같음. 그래도 두번째껀 실행하는 시간은 빨라서 망정이지 얘는... 저 모듈 Jupyter에서 불러오지도 못함 OTL
전체 코드
import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
# gedit으로 여는데 엄청 방대해서 랙걸린다... (리눅스긴 한데 PC 사양이 그렇게 좋은 편은 아님)
Chromosome=[]
CLNSIG=[]
Gene=[]
# 일단 Series를 거쳐서 DataFrame화 할 예정.
for record in vcf_reader:
INFO_get=record.INFO.get('CLNSIG')
Gene_get=record.INFO.get('GENEINFO')
if INFO_get:
Chromosome.append(record.CHROM)
CLNSIG.append(str(INFO_get))
else:
Chromosome.append(record.CHROM)
CLNSIG.append("No data")
# CLNSIG에 키가 없으면 처리하는거
if Gene_get:
Gene.append(Gene_get)
else:
Gene.append("No data")
# Gene도 몇개 없더라...
# 이거 자체는 간단하다. 키가 없으면 NaN으로 넣는다는 얘기.
Chr_series=pd.Series(Chromosome)
CLN_series=pd.Series(CLNSIG)
Gene_series=pd.Series(Gene)
# 제발 시리즈 한번에 나오게 해주소서
Gene_df=pd.DataFrame({"Chromosome":Chr_series, "Gene":Gene_series,"CLNSIG":CLN_series})
# 오케이 제발 데이터프레임 이쁘게 나오게 해주소서
Gene_df2=Gene_df.groupby(['Chromosome','CLNSIG']).count()
# Groupby로 시원하게 묶어보자
# 근데 이거 object라 정렬이 1, 2, 3, 4로 안된다...
Gene_df3=Gene_df.groupby('Gene').count()
Gene_df3=Gene_df3.sort_values('CLNSIG',ascending=0)
print(Gene_df3.head(30))
# 정렬하고 30개 뽑아보자
문제
Q1. VCF파일 내 데이터 중 CLNSIC 데이터를 염색체별로 분류할 것
Q2. Pathogenic한 Gene Top 30 꼽기
공통 절차-분석을 위한 Dataframe 만들기
import vcf
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
for record in vcf_reader:
print(record)
# PyVCF로 불렀는데 상당히 방대하다. 와... 이거 판다스 있어야된다...일단 저걸 하려면 VCF 파일을 열어야 한다. VCF file의 경우 ClinVar에서 받을 수 있고, ftp에 일주일마다 업그레이드가 된다. 아무튼... 저거는 바이오파이썬으로 불러오는 게 아니라 PyVCF라는 모듈이 또 있어서 설치를 해야 하는데 여기서부터 난관이었다.
1. 설치를 했는데 Jupyter에서 못 불러온다
2. 파이참에서도 못 불러온다
3. IDLE은 되네?
4. 어? 갑자기 파이참에서 되네? (Jupyter는 안됨)
참고로 코드 실행 속도가 정말 어마무시하게 느리다. 데이터가 방대하니 어쩔 수 없지...
import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
for record in vcf_reader:
print(record.INFO['CLNSIG'])
# PyVCF로 불렀는데 상당히 방대하다. 와... 이거 판다스 있어야된다...
CLNSIG_frame=pd.DataFrame(record.INFO['CLNSIG'])그리고 그냥 이렇게 Dataframe 만들려고 하면 Key error가 난다. 찾아보니 해당하는 키가 없으면 나는 에러라고 한다.
import vcf
import pandas as pd
# 모듈은 항상 위쪽에 부릅니다.
vcf_reader = vcf.Reader(open('/home/koreanraichu/clinvar_20211204.vcf', 'r'))
# gedit으로 여는데 엄청 방대해서 랙걸린다... (리눅스긴 한데 PC 사양이 그렇게 좋은 편은 아님)
Chromosome=[]
CLNSIG=[]
# 일단 Series를 거쳐서 DataFrame화 할 예정.
for record in vcf_reader:
INFO_get=record.INFO.get('CLNSIG')
if INFO_get:
Chromosome.append(record.CHROM)
CLNSIG.append(INFO_get)
print(record.CHROM,INFO_get)
else:
Chromosome.append(record.CHROM)
CLNSIG.append("NaN")
print(record.CHROM,"No keys")
# 이거 자체는 간단하다. CLNSIG에서 Key error가 나는데, 없으면 NaN으로 쓴다는 얘기.그래서 이렇게 if함수를 넣어줍니다. Key가 있으면 해당하는 값을 넣고, 없으면 Nan(현재는 No data)을 넣으라는 얘기.
Chr_series=pd.Series(Chromosome)
CLN_Series=pd.Series(CLNSIG)
# 심보러시여 제발 시리즈 한번에 나오게 해주소서에이 솔직히 시리즈는 잘 돼요... 면접관님 보시는데 무슨 심보러 드립을 치고 있냐
CLN_df=pd.DataFrame({"Chromosome":Chr_series, "CLNSIG":CLN_Series})
print(CLN_df)
# 오케이 제발 데이터프레임 이쁘게 나오게 해주소서데이터프레임...... (마른세수)
중간에 개고생하긴 했으나 해결.
Chromosome CLNSIG
0 1 [Uncertain_significance]
1 1 [Likely_benign]
2 1 [Benign]
3 1 [Likely_benign]
4 1 [Uncertain_significance]
... ... ...
1105836 MT [Uncertain_significance]
1105837 MT [Uncertain_significance]
1105838 MT [Benign]
1105839 MT [Uncertain_significance]
1105840 NW_009646201.1 [Benign]
[1105841 rows x 2 columns]
Process finished with exit code 0이렇게 데이터프레임이 나왔다. (참고로 최종 Dataframe은 2번 문제도 하나의 데이터프레임으로 풀기 위해 Gene 컬럼도 추가되어있음)
CLNSIG 집계
CLN_df2=CLN_df.groupby('Chromosome').count()
print(CLN_df2)CLNSIG
Chromosome
1 89545
10 41010
11 65665
12 49168
13 29113
14 33996
15 38488
16 60126
17 76248
18 18665
19 51709
2 112429
20 20230
21 12514
22 21575
3 56286
4 34468
5 58692
6 48407
7 51289
8 36992
9 50025
MT 2838
NW_009646201.1 1
X 46316
Y 46참고로 저게 Object라 정렬이 저렇게 된다... 예전에 scatter plot 그릴때도 저랬는데, 그때는 결측값 지우고 object에서 numeric으로 바꿔서 해결 봤다... 그건 축 값이 순수하게 숫자라 그런거고 저거는 염색체라 지우는 게 안돼요... (주륵)
CLN_df2=CLN_df.groupby('Chromosome').aggregate(['CLNSIG'])
print(CLN_df2)
# 염색체별로 묶어보자.이거 뭐 계속 바꿔서 해봤더니 아 리스트라 안된다고!!! 하더라... 그 답이 전체 코드에 있다.
CLNSIG.append(str(INFO_get))정확히는 이거. 저 INFO에 CLNSIG 말고 데이터가 많은데 그걸 다 리스트로 불러와서 저런 사태가 일어난 것이라고 한다. 딕셔너리 키값은 튜플이나 문자열같이 절대불변인 것들로 해야 하기 때문. 그래서 if문 돌리고 리스트에 넣을 때 문자열로 변환했다.
Gene_df2=Gene_df.groupby(['CLNSIG','Chromosome']).count()
print(Gene_df2)
# Groupby(CLNSIG, Chromosome 순)Gene
CLNSIG Chromosome
No data 1 56
10 29
11 57
12 18
13 10
... ...
['risk_factor'] 6 30
7 16
8 21
9 20
X 10
[527 rows x 1 columns]니네 근데 진짜 이걸로 해결되는거 실화냐고... (No data: CLNSIG 데이터가 없음)
Gene_df2=Gene_df.groupby(['Chromosome','CLNSIG']).count()
print(Gene_df2)
# GroupbyGene
Chromosome CLNSIG
1 No data 56
['Affects'] 19
['Benign', '_other'] 1
['Benign'] 15397
['Benign/Likely_benign', '_other'] 2
... ...
Y ['Benign'] 4
['Likely_benign'] 5
['Likely_pathogenic'] 2
['Pathogenic'] 28
['Uncertain_significance'] 7
[527 rows x 1 columns]난 저거 말고 염색체별로 볼 건데? 하면 순서를 바꾸시면 됩니다.
Top 30 pathogenic gene
이거 Top 30 뽑는거는 간단하다. ascending에 0 주고 head(30)으로 뽑으면 된다.
Gene_df3=Gene_df.groupby('Gene').count()
Gene_df3=Gene_df3.sort_values('CLNSIG',ascending=0)
print(Gene_df3.head(30))
# 정렬하고 30개 뽑아보자Chromosome CLNSIG
Gene
BRCA2:675 13480 13480
BRCA1:672 11565 11565
TTN:7273|TTN-AS1:100506866 9999 9999
APC:324 8901 8901
NF1:4763 7658 7658
TTN:7273 7601 7601
TSC2:7249 6392 6392
ATM:472 6260 6260
MSH6:2956 5591 5591
POLE:5426 4939 4939
FBN1:2200 4796 4796
RYR2:6262 4358 4358
MSH2:4436 4262 4262
RYR1:6261 4063 4063
DMD:1756 4063 4063
ATM:472|C11orf65:160140 3829 3829
PALB2:79728 3788 3788
NEB:4703 3706 3706
SYNE1:23345 3486 3486
DICER1:23405 3403 3403
MLH1:4292 3354 3354
BRIP1:83990 3350 3350
USH2A:7399 3300 3300
PLEC:5339 3267 3267
SMARCA4:6597 3024 3024
PMS2:5395 2949 2949
LDLR:3949 2842 2842
TSC1:7248 2776 2776
POLD1:5424 2736 2736
CDH1:999 2715 2715이거는 유전자 전체 통틀어서 출력해주는거고
Gene_df3=Gene_df.groupby(['Gene','CLNSIG']).count()
Gene_df3=Gene_df3.sort_values('Chromosome',ascending=0)
print(Gene_df3.head(30))
# 그냥 Gene으로 정렬하는 코드는 이거Chromosome
Gene CLNSIG
BRCA2:675 ['Uncertain_significance'] 5467
APC:324 ['Uncertain_significance'] 4835
TTN:7273|TTN-AS1:100506866 ['Uncertain_significance'] 3933
BRCA2:675 ['Pathogenic'] 3590
TTN:7273 ['Uncertain_significance'] 3238
ATM:472 ['Uncertain_significance'] 3207
NF1:4763 ['Uncertain_significance'] 3101
POLE:5426 ['Uncertain_significance'] 3093
BRCA1:672 ['Uncertain_significance'] 3012
MSH6:2956 ['Uncertain_significance'] 2944
BRCA1:672 ['Pathogenic'] 2892
TTN:7273|TTN-AS1:100506866 ['Likely_benign'] 2680
BRCA1:672 ['not_provided'] 2654
BRCA2:675 ['Likely_benign'] 2445
TTN:7273 ['Likely_benign'] 2441
TSC2:7249 ['Uncertain_significance'] 2212
RYR2:6262 ['Uncertain_significance'] 2123
MSH2:4436 ['Uncertain_significance'] 2002
ATM:472|C11orf65:160140 ['Uncertain_significance'] 1995
PALB2:79728 ['Uncertain_significance'] 1956
BRIP1:83990 ['Uncertain_significance'] 1917
NF1:4763 ['Pathogenic'] 1895
RYR1:6261 ['Uncertain_significance'] 1757
DICER1:23405 ['Uncertain_significance'] 1745
SYNE1:23345 ['Uncertain_significance'] 1712
PLEC:5339 ['Uncertain_significance'] 1670
BRCA1:672 ['Likely_benign'] 1617
NEB:4703 ['Likely_benign'] 1617
NF1:4763 ['Likely_benign'] 1596
RECQL4:9401 ['Uncertain_significance'] 1577
Process finished with exit code 0이렇게 하면 유전자별로 가장 높은 CLNSIG만 나온다.
'Coding > Python' 카테고리의 다른 글
| 10진수->2진수 변환 코드 (0) | 2022.08.21 |
|---|---|
| 번외편-코딩테스트 풀이 (3) (0) | 2022.08.21 |
| Biopython-dbSNP와 Clinvar (0) | 2022.08.21 |
| 심심해서 써보는 본인 개발환경 (0) | 2022.08.21 |
| Biopython-Q&A (0) | 2022.08.21 |