2108 Hakyung TUXEDO

From Crop Genomics Lab.
Jump to: navigation, search

Contents

Index Building

tophat2는 fasta file을 그대로 사용하기 보다는 index를 만들어 파일 사용을 용이하게 한다. 이를 위해 첫 턱시도 파이프라인인 tophat을 사용하기 앞서 bowtie2-build index를 사용한다. command 는 다음과 같으며, output directory name은 웬만하면 referece.fa 와 같은 이름으로 해주는 것이 나중에 헷갈리지 않고 좋다.

bowtie2에 대한 자세한 설명은 링크 참조.

bowtie2-build <reference.fa> <output dir name>

Tophat

tophat은 read들을 mapping 하는 과정이다. 각 sample들을 각각 mapping해주며, pair-end 의 경우 아래와 같은 command를 입력해주면 된다. output은 -o 뒤에 지정해준 이름으로 폴더 안에 저장되며, accepted.bam이 앞으로 사용해야할 bam file이다. ref_index에는 위에 bowtie2-build index에서 지정해주었던 이름을 적으면 된다. -p 는 사용할 thread의 갯수인데, htop을 사용하여 서버 상태를 보고 정할 것. bowtie2 index를 사용했을 경우, tophat2만 되던데... 244에서만 그런건지 원래 그런건지 잘 모르겠음. ㅠ

tophat2 -p 2 -o <*> -G <anotation.gff3> <ref_index> <*_1.fastq.gz> <*_2.fastq.gz>

cufflinks

cufflinks -o <*_cufflinks> -G <annotation.gff3> <*/accepted_hits.bam>

HTseq-count

alined 된 read들을 세어준다. cufflink를 이용한 TUXEDO pipeline을 사용하는 경우, 1:1 비교를 하므로 replication이 많을 경우에는 모든 경우의 수를 생각해주거나, 아니면 반복적으로 실행한 후에 일정 조건을 충족시키는 애들만 남겨놓거나... 아무튼... 복잡하므로 HTseq을 이용해 주는 것이 좋다. HTseq이 요구하는 pipeline을 처음부터 따르지 않은 경우에는 --stranded=no 옵션을 반드시 넣어주어야 한다. 기본 input이 sam이므로 만일 input file에 bam이 아닌 sam을 넣어주고자 한다면 -f 옵션은 넣어주지 않아도 된다. annotation gtf file은 기존에 사용한 gff3 file과 내용은 동일하나 포맷만 바꿔준 것.

자세한 사항은 HTseq 링크를 참조.

python -m HTSeq.scripts.count -t exon -f bam stranded=no <*_sorted.bam> <annotation.gtf>  > *_HTSeq.txt

edgeR

edgeR을 실행하기 앞서 HTseq에서 세어준 결과들을 한 테이블에 모아 주어야 한다. 각 coumn들이 sample들의 count수를 포함하게 병렬로 병합하면 된다. 결과 중에 logFC (log2 fold change) 값이 1보다 크거나 (2배 이상 증가), -1 보다 작은 것 (2배 이상 감소), 그리고 P value가 0.05 이하인 것들을 보통 DEG로 간주한다.

Library(edgeR)
x<- read.delim(“m1.m2.f1.f2.htseq.count(합친파일)”,row.names=“Symbol-”) : or 미리 table에 tag를 달아주었다면 read.table(변수)로 처리해도 됨.
group <- factor(c("M","M","F","F"))
y <- DGEList(counts=x,group=group)
y <- calcNormFactors(y)
design <- model.matrix(~group) (If you have two trait, ~trait1+trait2)
y <- estimateDisp(y,design)
fit <- glmQLFit(y,design) 
qlf <- glmQLFTest(fit) (통계적으로 유의미한 것 고르기)
topTags(qlf)
write.table(x=qlf$table, file="m1.m2.f1.f2.htseq.count.edgeR") output 저장