cnvkit 설치 방법

간만에 프로그램 설치 기록 하나 남겨 봅니다.

설치할 tool은 Copy Number Variant 탐지에 대명사 중 하나인 CNVKit입니다.

참조할 사이트는 당연히 공식 cnvkit github고요 

첫번째 방법은 conda를 이용한 방법!

$ conda config --add channels defaults 
(conda를 설치해서 사용하시는 분들은 굳이 defaults는 안하셔도 됩니다.)

$ conda config --add channels conda-forge
$ conda config --add channels bioconda
$ conda install psutil
$ conda create -n cnvkit cnvkit
라고 명령어를 입력하면 엄청 많은 패키지들을 설치하겠다고 나옵니다.

$ source activate cnvkit

(cnvkit)$ python cnvkit.py -h

numpy가 없다고 나와서 다음과 같이 가상 환경 안에서 numpy 설치

(cnvkit)$ conda install numpy

그러나 conda를 이용한 방법은 제대로 작동을 안해서...

$ source deactivate

이유는 꼭 알아야 할까요? 다른 방법으로 설치하면되지요!

그래서 걍 github에서 소스 파일 다운 받아서 pip를 사용하여 설치!!

$ git clone https://github.com/etal/cnvkit

$ cd cnvkit/

$ pip install -e .

$ ~/Python-2.7.12/bin/python cnvkit.py -h

기본 pip를 사용했기 때문에 python-2.7.12에 설치가 되었고, 만약 pip3를 사용하여 설치하였다면... python3을 이용하면 되었을 것으로 예상됩니다.

그래서 오늘은 conda와 pip를 이용한 설치 방법에 대해서 알아보았습니다.

다음에 또 유익한 정보를 가지고 찾아오도록 하겠습니다.

출처: @candyz_hyojung

NCBI기반의 분류데이터를 손쉽게 다루는 방법

NCBI의 taxonomy 데이터는 드려다 보신 분들은 아시겠지만 환장의 도가니라고 할 수 있죠

그래서 보다 못한 호주의 어느 용자께서 NCBI의 taxonomy 데이터를 처리할 수 있는 스크립트를 개발하셔서 github에 올려주셨습니다. 물론 논문도 투고하셨고요

이름하여 ncbi-taxonomist, NCBI 분류학자라니.. Orz

Collecting and managing taxonomic data with NCBI-taxonomist

doi: https://doi.org/10.1093%2Fbioinformatics%2Fbtaa1027

여튼 재미있는 기능 중 하나가 NCBI의 nucleotide이든 protein이든 accession을 입력하면 이 accession의 taxonomy 정보를 알려준다는 것입니다.

물론 NCBI 분류학자답게 NCBI의 taxonomy ID과 종 이름(약간의 오타가 있어도, 그 오타가 NCBI에 등록되어 있다면)로도 검색할 수 있습죠.

accession 정보를 이용하여 taxonomy 정보를 파악하는게 생각보다 번거로운 작업인데 ncbi-taxonomist에서 명령어 한번 때려주면 호로록 검색해서 결과를 알려줍니다.

여튼 오늘은 간단한 기능을 보여주는 tool을 소개하였는데
밤새 태풍에 무탈하시기 바랍니다. :)

출처: @ye._.vely618

AAPL 1년치 주식 자료를 쉽게 다운로드 받을 수 없을까?

왜 없겠습니까?

그래서 chatGPT와 Bing에게 물어보았습니다.

chatGTP

----

Q: https://finance.yahoo.com/ 에서 API를 이용해서 AAPL 주식 1년 자료 다운로드 받는 코드 짜줘.

A: 안녕하세요! Yahoo Finance API는 2017년 이후에 deprecated되어 더 이상 지원되지 않기 때문에, 이를 이용한 데이터 수집이 불가능합니다. 대신 Yahoo Finance에서는 다양한 주식 데이터를 제공하고 있으므로 다른 방법으로 데이터를 수집해야 합니다.

먼저, Python에서는 yfinance 라이브러리를 이용해서 Yahoo Finance에서 주식 데이터를 다운로드할 수 있습니다. yfinance를 사용하면 간편하게 주식 데이터를 다운로드할 수 있습니다. 아래는 yfinance 라이브러리를 이용해 AAPL의 1년치 주식 데이터를 다운로드하는 코드입니다.

import yfinance as yf

import datetime

start_date = datetime.datetime.now() - datetime.timedelta(days=365)

end_date = datetime.datetime.now()

# AAPL의 주식 데이터를 다운로드합니다.

aapl = yf.download("AAPL", start=start_date, end=end_date)

# CSV 파일로 저장합니다.

aapl.to_csv("AAPL.csv")

이 코드는 현재 날짜를 기준으로 1년 전부터 현재까지의 AAPL 주식 데이터를 다운로드하고, CSV 파일로 저장합니다. 위 코드에서 start_date와 end_date 변수를 변경함으로써 원하는 기간의 데이터를 다운로드할 수 있습니다.

----

BING

----

# yfinance 라이브러리 설치

pip install yfinance

# yfinance 라이브러리 임포트

import yfinance as yf

# AAPL 티커 객체 생성

aapl = yf.Ticker("AAPL")

# 1년 동안의 주가 데이터 가져오기

hist = aapl.history(period="1y")

# 데이터를 csv 파일로 저장하기

hist.to_csv("aapl_1y.csv")

----

chatGTP와 bing에게 코드를 작성해 달라고 하였을 때 현재 단순한 코드에서는 chatGPT의 손을 들어주고 싶습니다.

그 이유는 date와 period의 차이라고 해두죠. 나중에 특정 년 월 일을 사용할 것을 생각한다면 bing보다는 chatGPT가 (아직은) 한 수 위 인듯 합니다. :)

출처: @ye._.vely618

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

Biopython Entrez 사용하기

최근에 rosalind를 다시 해보면서
가능하면 파이썬 기본 라이브러리로 해결해보려고 했는데
Bioinformatics Armory 부분의 일부 문제는 Biopython을 사용하지않으면
꽤나 곤란한 상황이 발생할것 같아서
일부 문제에서는 걍 추천하는대로 biopython을 사용하는걸로..

그리고 Biopython과 함께 E-utilities를 사용해야 하는데
역시 간만에 보니 생각이 안나는게 인지상정

관련 NCBI 페이지


E-utility에서 사용하는 DB 이름..

Entrez Database	UID common name	E-utility Database Name
BioProject	BioProject ID	bioproject
BioSample	BioSample ID	biosample
Biosystems	BSID	biosystems
Books	Book ID	books
Conserved Domains	PSSM-ID	cdd
dbGaP	dbGaP ID	gap
dbVar	dbVar ID	dbvar
Epigenomics	Epigenomics ID	epigenomics
EST	GI number	nucest
Gene	Gene ID	gene
Genome	Genome ID	genome
GEO Datasets	GDS ID	gds
GEO Profiles	GEO ID	geoprofiles
GSS	GI number	nucgss
HomoloGene	HomoloGene ID	homologene
MeSH	MeSH ID	mesh
NCBI C++ Toolkit	Toolkit ID	toolkit
NCBI Web Site	Web Site ID	ncbisearch
NLM Catalog	NLM Catalog ID	nlmcatalog
Nucleotide	GI number	nuccore
OMIA	OMIA ID	omia
PopSet	PopSet ID	popset
Probe	Probe ID	probe
Protein	GI number	protein
Protein Clusters	Protein Cluster ID	proteinclusters
PubChem BioAssay	AID	pcassay
PubChem Compound	CID	pccompound
PubChem Substance	SID	pcsubstance
PubMed	PMID	pubmed
PubMed Central	PMCID	pmc
SNP	rs number	snp
SRA	SRA ID	sra
Structure	MMDB-ID	structure
Taxonomy	TaxID	taxonomy
UniGene	UniGene Cluster ID	unigene
UniSTS	STS ID	unists

모 그렇다고 합니다.

sensors 로그 기록용 스크립트

서버 온도가 아무 작업을 안해도
60도를 상회하는 관계로 sensors로 온도 변화를 기록하고자
급하게 작성한 한 스크립트.
물론 쓸모는 없었다는게 함정 ㅋ

import time
from datetime import datetime
from subprocess import check_output

while 1:
        print datetime.today().strftime("%Y/%m/%d/ %H:%M:%S")
        contents = check_output(['sensors'])
        print contents
        time.sleep(60)

DELL 서버들은 iDRAC에서 로그를 뽑아내서
엔지니어한테 전달해줘야봐준다능... Orz

파일의 포맷을 변환하는데 필요한 것들

내가 아니란 말이닷!!! ㅋㅋ

python에서 Biopython을 이용하여
간단하게 convert하는 샘플 코드를 제공하고 있으니
여러분들도 쉽게 만들수 있어요~ :)
Biopython에서 제공하는 Tutorial 


오늘 문의가 들어온 파일은 sff파일
Roche의 454 GS FLX? sequencing 결과파일로....
ABI와 함께 illumina한테 밀려서 뒷방으로 들어앉은 파일 포맷입니다.
그러나 아직도 쓰는 이유는 read 길이가 길기때문 :)

그렇습니다. PacBio도 Nanopore다 디립다 길게 sequencing해준다는
애들이 있습니다. 그런데 왜 옛날꺼 쓰냐?? PacBio는 base quality가 안습이고,
Nanopore는.... 언제 출시일지 전 잘 모르겠습니다. 업자가 아닌관계로 ㅎㅎ

그래서 위의 길게 sequencing 해준다는 시퀀서를 제외하고는 Roche의 454가 read 길이가 가장 길다고 할 수 있겠습니다. NGS중에선 말이죠

그런데 sff파일을 보려고 하면 문제가 생깁니다.
권모씨께서 문의를 한것이 그것때문인지는 모르겠지만 걍 일반인이
sff파일을 걍 직접 볼수가 없습니다. 왜냐 binary파일이니깐요(sff파일이 binary라고
알고 있는데  직접 다뤄본적이 없어서... ㅎㅎ )

그래서 사람이 볼수 있게 파일을 변환시켜줘야 한다는 겁니다.

convertSff.py

#!/usr/bin/python

import os, sys
from Bio import SeqIO

try:
inputSFF = sys.argv[1]
outputPREFIX = sys.argv[2]

except:
print "Usage: python convertSFF <input.sff> <output_name>"
print ""
exit(1)


SeqIO.convert(inputSFF,"sff","%s.fasta"%(outputPREFIX), "fasta")
SeqIO.convert(inputSFF,"sff","%s.quality"%(outputPREFIX), "qual")
SeqIO.convert(inputSFF,"sff","%s.fastq"%(outputPREFIX), "fastq")

권모씨의 요청으로 급조한 날림 convert python 코드 ㅋㅋ
이 스크립트를 수행하면 세개의 파일이 나오게 될것으로 예상됩니다. ㅎㅎ
안나오면 어쩔수없고... ㅎㅎ


아.. 그리고 사족으로 LT사의 SOLiD의 경우 우리가 알고 있는 서열과 달리
첫 염기 서열만 서열이고 그 다음부터는 A/G/T/C 알파벳이 아닌 숫자로 되어있는데..
이걸 굳이 변환해서 reference geneome에 mapped 작업하지 말라고 합니다.
Re-sequencing하는 경우라면 변환해서 mapping하지 말고 원래 원본 파일 그대로를
input으로 하는 align 프로그램을 사용해서 mapped한 다음에 그 다음 작업을
일반적으로 사용하는 samtools나 GATK같은 프로그램을 사용하라고 합니다.
(다들 알고있는거 한번더 상기 시켜드렸습니다. 혹시 아나요 SOLiD 포맷을 분석하게 될지.. ㅎㅎ)

분석시 raw 파일을 사용해야 하는 이유는 SOLiD만의 월등한 quality 효과를 볼수 있어서
그러지 않겠나하는....  믿거나 말거나 저 혼자만의 생각입니다.. ㅎㅎ
다만, 타사 제품과 다르게 복잡하게 숫자로 표현한건 아니겠죠...
나름의 숨은 뜻이.... 쿨럭.. (설마... 간지용;;;;; )

Re-sequencing이 아닌 denovo일 경우 모 어쩔수 없이 fastq파일로 변환을 해야 하지 않을까 합니다. assembly 프로그램을 작동시키려면 아무래도 SOLiD format보다는 fastq 포맷이
수월하니깐요.. :)

그럼....

Python에서 통계 작업하기

통계작업할 때 지금까지는 대부분 엑셀,
조금 데이터양이 많아지거나 고급 통계기법이 필요하면 R을 사용하였다.

근데 아쉽게도 R에서의 문자열 핸들링관련해서 지식이 민망하여
통계처리 중간에 문자열 조작이 필요한 경우 중간에 값을 텍스트로 빼서
파이썬에서 작업하고 다시 R로 가져가는 상당히 귀찮은 작업을 거쳤다.

그래서..
결국 미루다 미루다..

통계작업 중간에 문자열작업이나 다른 수정이 필요업는 단계까지
파이썬에서 작업을 마치는 방법을 터득하기로 했다.

근데.. 이게 파이썬에서 통계작업을 한번도 안해본지라
몇날 몇일 헤매고 있었는데

Scipy(파이썬 설치할때 기본적으로 설치하는 모듈인데... 쓰지는 않고 있었다)에서
통계 함수들을 제공한다니.. ;;;

Scipy Statistical Function

여하튼...
이것때문에.. 어느정도 고민 해결~ ㅎㅎ

Fasta 형식의 파일에서 빈서열 제거

가끔씩 Blast를 수행하고자 formatdb를 수행 할 때,
다음과 같은 에러를 접한 적이 있으리라 본다.

[formatdb] WARNING: Cannot add sequence number XXXXX XX.XXX.XXX.
 because it has zero-length.
[formatdb] FATAL ERROR: Fatal error when adding sequence to BLAST database.

formatdb를 수행하려는 fasta 서열에
빈 서열을 가지고 있기 때문에 나오는 에러로
빈 서열을 제거하면 OK!

vi check.py

import glob,sys
from Bio import SeqIO


file = sys.argv[1]


ow = open(file.split('.')[0]+'.check','w')
seqs = SeqIO.parse(open(file), format='fasta')
 for rec in seqs:
        name = rec.description
        seq = rec.seq.tostring()


        if len(seq.strip()) != 0:
                ow.write('>'+name+'\n')
                ow.write(seq+'\n')


ow.close()

[test]$python check.py test.fasta

Win32com Python

파이썬에서 좀더 다양한 윈도우 프로그래밍이 가능하게 하는 모듈?
여하튼.... 좋은거~ ㅋ

Win32com

다양한 윈도우 프로그램과 연동되어
프로그래밍을 더 시키게 만든 녀석.. ㅡ.-a

여하튼.... python 2.3이후부터는 없는줄 알았는데..
역시 sourceForge에는 있었다는...;;;

대용량 파일을 위한 shelve

출처: 파이썬마을


파이썬에 이런 모듈이..
참.. 대단한듯...
-파이썬에 감탄하면 지는거다..
난 자바를 공부하고 있어야돼... ㅠ.ㅠ

쉽게말해 대용량 사전
-사전이라는 자료형 은근히 편했는데..


실제 사용 예로 EMBL CDS fasta파일을 읽어서 E. coli에서 나온 것만 골라서
뽑아주는 소스는 이렇게 만들 수 있습니다.

from Bio import SeqIO
import shelve

seqs = SeqIO.parse(open('/data/embl/cds_nr.fasta'), format='fasta')
db = shelve.open('ecoli.db')


for rec in seqs:
    if ('Escherichia coli' in rec.description and
        'hypothetical' not in rec.description):
            name = rec.description.split()[0].split(':')[1]
            seq = rec.seq.tostring()
            db[name] = seq

나중에 EMBL ID로 검색할 때는 다음과 같이 open만 하면 검색할 수 있습니다.

import shelve

db = shelve.open('ecoli.db')
genes = ['100746833', '1008171518']
for gene in genes:
    print db[gene]

32G가 30여분 소요 됐다니.. 난 4.1G(nr, 한국시간 2009.3.5에 다운로드) 밖에 안돼네.. ㅋㅋ

5분이면 끝나야 하는거 아니야..;;
왜 안끝나..;;

Python 64Bit Library

예전에는 라이브러리를 사용하는데 있어서
요구되는 라이브러리들의 버전만 체크하면 되는데
요즘 64bit OS들이 널리 퍼지면서(본인 또한 64bit OS만 쓴다는.. 쿨럭...)
버전 체크 뿐만 아니라 bit확인도 해주어야 한다.

그런데 python을 사용하다 보면 32bit 버전 라이브러리는 홈페이지에서
다운받기 쉬운데 64bit 버전 라이브러리는 좀처럼 찾기 쉽지 않은 것들이 있다.

google 검색하다보니 다음과 같이 64bit용 python 라이브러리를
모아둔곳이 까먹기 전에 게시~ :)

http://www.lfd.uci.edu/~gohlke/pythonlibs/

Python에서 집합 연산

파이썬에서 집한연산을 하는 경우가 생기게된다.

이럴때 쓰라고 있는 set 

집합형으로 바꾸는 방법은 걍 set()을 묶어주면 Ok.

list > set

리스트의 원소 중 중복되는 원소는 제거하고 집합형으로 변환

string > set

문자열도 한단어씩 split후 중복되는 단어 제거하고 집합형으로 변환

교집합

setA.intersection(setB)

set.intersection(setA,setB,setC)

합집합

setA.union(setB)

set.union(setA,setB,setC)

차집합

setA.difference(setB)

set.difference(setA,setB,setC)

해석: setA.difference(setB.union(setC))

대칭차집합 (두 set에 대해서만 계산 가능)

setA.symmetric_difference(setB)

set.symmetric_difference(setA,setB)

setA = set([1,2,3,4,5,6])

setB = set([1,2,3])

부분집합(subset)

옳은예: setB.issubset(setA) ( = set.issubset(setB,setA)

틀린예: setA.issubset(setB) ( = set.issubset(setA,setB)

앞쪽 집합이 뒤쪽에 나오는 집합의 부분집합이냐고 확인

초집합(superset)

부분집합과 반대로 생각하면됨

옳은예: setA.issubset(setB) ( = set.issubset(setA,setB)

틀린예: setB.issubset(setA) ( = set.issubset(setB,setA)