반응형

Argument받아서 처리하기




bash 스크립트로 argument를 여러 개 받아 해당 수 만큼 for loop을 돌릴 때 쓰는 방법이다.


전체 argument를 반복하고 싶다면 아래처럼 하면 된다.


for i in $@ ; do echo $i ; done 


특정 구간만을 반복시키고 싶다면 indexing을 넣어주어야 한다. 아래 예시는 2번부터 끝까지 반복이다.


for i in ${@:2} ; do echo $i ; done


2번부터 3번 반복하고 싶다면 (즉 2부터 4까지) 아래처럼 더 추가한다.


for i in ${@:2:3} ; do echo $i ; done



argument가 몇 개 인지 알고싶다면 "$#" 를 사용하면 된다.


argument가 0개 일때의 조건문.


if [ $# -eq 0 ] ; then

        echo "Number of argument is 0"

else 

        echo "Number of argument is $#"

fi





Reference -

https://stackoverflow.com/questions/255898/how-to-iterate-over-arguments-in-a-bash-script

https://stackoverflow.com/questions/4423306/how-do-i-find-the-number-of-arguments-passed-to-a-bash-script

반응형
반응형

특정 파일이나 디렉토리를 남기고 모두 지우기



file.txt라는 이름의 파일을 남기고 나머지를 모두 지우기.


find . ! -name 'file.txt' -type f -exec rm -f {} +



디렉토리라면 type을 d로 바꾸고 rm 에 recursive 옵션을 추가하면 된다.


find . ! -name 'file.txt' -type d -exec rm -rf {} +

* 주의사항 

위의 디렉토리만 남기고 모두 지우기는 디렉토리 안의 파일도 모두 지움

디렉토리 수준에서만 적용하고 싶다면 아래처럼 사용한다


find . -maxdepth 1 ! -name 'file.txt' -type d -exec rm -rf {} +
depth제한을 둠으로써 하위 폴더로 내려가지 않고 현재 폴더에서만 적용.


여러 개의 파일을 남겨놓고 싶다면 아래처럼 추가 할 수 있다.


find . ! -name 'file.txt' ! -name 'file2.txt' -type f -exec rm -f {} +


Reference -

https://unix.stackexchange.com/questions/153862/remove-all-files-directories-except-for-one-file

반응형

'Computer Science > linux' 카테고리의 다른 글

linux 계정 정보 옮기기  (1) 2019.08.21
tar 디렉토리 지정해서 압축 풀기  (0) 2018.10.04
리눅스에서 프록시 설정하기  (0) 2018.07.26
neocomplcache vim plugin 설치하기  (0) 2018.07.12
samba 설정하기  (0) 2018.02.02
반응형

16진수 RGB코드 알아내는법




화면 내에서 특정 색이 가진 RGB코드를 쉽게 알아내려면 EST soft에서 제공하는 알씨 프로그램이 가장 편하다. (https://www.altools.co.kr/download/alsee.aspx)


알씨 프로그램에서 [보기]->[픽셀정보] 를 클릭하면 마우스의 위치를 따라가면서 RGB값을 표시해준다. 0x는 16진법이라는 표시니 제외하고 뒤에 숫자+알파벳 조합이 색의 정보이다.



반응형
반응형

DESeq2에서 heatmap, PCA, MA, volcano plot 그리기




분석을 위한 DataSet이 필요하다. 예시에서는 여러 cell line에 몇 가지 drug를 각기 다른 dose로 처리했을 때의 결과를 분석하고자 한다.


DEG분석은 dose에만 초점이 맞춰지도록 디자인 하였다.


res변수는 결과가 adjust p-value 낮은순서대로 sorting하여 저장하였다.


sampleTable<-data.frame(sampleName=sampleFiles, fileName=sampleFiles, dose=I(dose), drug=drug, cell_line=cell_line)

ddsHTSeq<-DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory=directory, design=~dose)

dds<-DESeq(ddsHTSeq)

res <- results(dds)

res <- res[order(res$padj),]


heatmap을 그리기 위한 color를 지정해 준다. 지정하지 않으면 그림을 그릴 때 마다 다른 color가 나오기 때문에 결과를 정리하는데 매 번 수정하기가 상당히 귀찮기 때문이다. 각 항목에 맞는 값을 넣어주고 뒤에 16진법 RGB코드를 넣어주면 된다. 

16진법 RGB코드는 아래 포스팅한 방법으로 찾으면 쉽게 찾을 수 있다.

## annotation color

ann_colors = list(

        dose = c(".." = "#F9EBD1"),

        cell_line = c(".." = "#FFFFFF"),

        drug = c(".." = "#BDBE32")

)


heatmap을 그리기 위해서 pheatmap 라이브러리를 로딩한다. order뒤에 [1:1000]이라고 되어 있는 부분은 유전자가 너무 많아서 발현량 기준으로 상위 1000개만 가져온 것이다. 전부 넣으려다가는 메모리 덤프가 생길 것이다. 1000개라는 갯수는 임의로 설정한 것이다.


## heatmap

library("pheatmap")

ntd <- normTransform(dds)

select <- order(rowMeans(counts(dds,normalized=TRUE)),

                decreasing=TRUE)[1:1000]

df <- as.data.frame(colData(dds)[,c("cell_line","drug","dose")])

rownames(df) = sampleFiles

pdf("heatmap_normalized_readcount.pdf")

pheatmap(assay(ntd)[select,], cluster_rows=T, show_rownames=F,

         cluster_cols=T, annotation_col=df, annotation_colors = ann_colors)

dev.off()


plotPCA는 함수 자체를 다시 정의한 것인데 이유는 DESeq2에서 제공하는 plotPCA함수는 PC1과 PC2밖에 고려하지 않기 때문이다.
마지막 ggplot을 불러오는 함수에서 x와 y를 PC3, PC4로 바꾸면 결과가 달라진다. PCA 분석에서 항상 PC1과 PC2의 조합이 정답이 아니기 때문에 필요한 방법이다.

## PCA
plotPCA <- function (object, intgroup = "condition", ntop = 500, returnData = FALSE)
{
    rv <- rowVars(assay(object))
    select <- order(rv, decreasing = TRUE)[seq_len(min(ntop,
        length(rv)))]
    pca <- prcomp(t(assay(object)[select, ]))
    percentVar <- pca$sdev^2/sum(pca$sdev^2)
    if (!all(intgroup %in% names(colData(object)))) {
        stop("the argument 'intgroup' should specify columns of colData(dds)")
    }
    intgroup.df <- as.data.frame(colData(object)[, intgroup,
        drop = FALSE])
    group <- if (length(intgroup) > 1) {
        factor(apply(intgroup.df, 1, paste, collapse = " : "))
    }
    else {
        colData(object)[[intgroup]]
    }
    d <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], PC3 = pca$x[, 3], PC4 = pca$x[, 4], group = group,
        intgroup.df, name = colnames(object))
    if (returnData) {
        attr(d, "percentVar") <- percentVar[1:2]
        return(d)
    }
    ggplot(data = d, aes_string(x = "PC1", y = "PC2", color = "group")) +
        geom_point(size = 3) + xlab(paste0("PC1: ", round(percentVar[1] *
        100), "% variance")) + ylab(paste0("PC2: ", round(percentVar[2] *
        100), "% variance")) + coord_fixed()
}
se <- SummarizedExperiment(log2(counts(dds, normalized=TRUE) + 1), colData=colData(dds,levels=c('0','0.1','1')))
pdf("PCA_normlaized_readcount.pdf")
plotPCA(DESeqTransform(se),intgroup=c("cell_line","drug","dose"),ntop=500)
dev.off()

MAplot은 유전자의 발현 값과 p-value의 상관관계를 보기 위해서 사용된다. 

#MAplot
pdf("MA_plot.pdf")
plotMA(res, ylim=c(-2,2))
dev.off()

Volcano plot은 각 유전자의 log2FoldChange와 p-value를 보기 위해 사용된다. plotting을 하는건 위의 4번째 줄 까지만 입력하면 되지만 특정 조건에 맞는 유전자일 때 색을 다르게 주거나 cutoff을 선으로 표시하기 위한 줄이 추가되어 있다.

#volcano plot
pdf("volcano_plot.pdf")

par(mar=c(5,5,5,5), cex=1.0, cex.main=1.4, cex.axis=1.4, cex.lab=1.4)
topT <- as.data.frame(res)
with(topT, plot(log2FoldChange, -log10(padj), pch=20, main="Volcano plot", cex=1.0, xlab=bquote(~Log[2]~fold~change), ylab=bquote(~-log[10]~Q~value)))

with(subset(topT, padj<0.05 & abs(log2FoldChange)>2), points(log2FoldChange, -log10(padj), pch=20, col="red", cex=0.5))

#Add lines for absolute FC>2 and P-value cut-off at FDR Q<0.05
abline(v=0, col="black", lty=3, lwd=1.0)
abline(v=-2, col="black", lty=4, lwd=2.0)
abline(v=2, col="black", lty=4, lwd=2.0)
abline(h=-log10(max(topT$pvalue[topT$padj<0.05], na.rm=TRUE)), col="black", lty=4, lwd=2.0)

dev.off()




반응형

'bioinformatics' 카테고리의 다른 글

oncotator 설치 및 실행하기  (0) 2018.10.04
liftover하기  (0) 2018.09.28
Optical duplicate와 Library duplicates  (0) 2018.08.27
Genome에 존재하는 Variant란?  (0) 2018.08.20
DESeq2에서 연속적인 값을 condtion으로 받기  (0) 2018.08.10
반응형

Cancer cell line 정보 받기




https://cancer.sanger.ac.uk/cell_lines


Catalogue Of Somatic Mutations In Cancer (COSMIC) 홈페이지에서 cancer cell line에 대한 정보를 받는 법을 알아보려고 한다. 특히나 특정 cancer cell line으로 실험을 하였을 때 결과가 제대로 나왔는지 확인하기 위해 해당 cell line에 존재하는 variant를 검사해 볼 수 있다.




가장 먼저 홈페이지에서 genome version을 설정해주어야 한다. 이 후의 자료들은 다 위의 genome version을 기반으로 제공된다.


옆에 검색창에 "HeLa"를 검색해보자.


비슷한 이름의 다른 정보가 하나도 없기 때문에 Samples에 하나 나온다. 나머지 정보는 뻔한 것들이니 넘어가고 뒷부분에 존재하는 variants를 보면 HeLa가 가지고 있는 variants가 보인다.



Filters를 사용하면 더 specific한 변이들만 따로 고를 수 있다.


Reference - 

https://cancer.sanger.ac.uk/cell_lines


반응형

'bioinformatics > cancer genomics' 카테고리의 다른 글

somatic mutation과 germline mutation  (0) 2018.10.12
Clinical Cancer 데이터베이스  (0) 2018.07.09
Molecular disease  (0) 2018.07.05
암 분류법  (0) 2018.07.05
CancerSCAN  (0) 2018.07.04

+ Recent posts