生物信息學編程實戰練習題大全

01-29

首先上目錄：

生信編程很簡單?

人類基因組的外顯子區域的長度?

hg19基因組序列的一些探究?

hg38每條染色體的基因、轉錄本分布?

多個同樣行列式文件的合併?

根據GTF畫基因的多個轉錄本結構?

下載最新版的KEGG信息，並且解析好?

寫超幾何分布檢驗?

ID轉換?

根據指定染色體及坐標得到序列?

根據指定染色體及坐標得到位置信息?

把文件內容按照染色體分開寫出?

JSON格式數據的格式化?

多個探針對應一個基因，取平均值、最大值?

把counts矩陣轉換成RPKM矩陣?

對有臨床信息的表達矩陣批量做生存分析?

對多個差異分析結果直接取交集並集?

根據GTF格式的基因注釋文件得到人所有基因的染色體坐標?

列出我們板塊的人氣最旺的20個題目，作為第五章節學習後的複習題目。

生信編程很簡單^1

編程語言系統入門

生信分析人員如何系統入門python?
生信分析人員如何系統入門perl？
生信分析人員如何系統入門R？
生信分析人員如何系統入門Linux？

題目

對FASTQ的操作

5,3段截掉幾個鹼基
序列長度分布統計
FASTQ 轉換成 FASTA
統計鹼基個數及GC%
對FASTA的操作
取互補序列
取反向序列
DNA to RNA
大小寫字母形式輸出
每行指定長度輸出序列
按照序列長度/名字排序
提取指定ID的序列
隨機抽取序列
高級難度
根據坐標取序列
多文件合併
根據ID列表取序列
GTF文件探索
簡併鹼基的引物序列還原成多條序列
snp進行注釋並格式化輸出

下載安裝BOWTIE2（內含測試數據）

cd ~/biosoftnmkdir bowtie && cd bowtienwget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip nunzip bowtie2-2.2.9-linux-x86_64.zipn

人類基因組的外顯子區域的長度^2

題目

下載人類外顯子的坐標文件，編寫代碼統計外顯子區域的長度。

測試數據

Rbioconductor的TxDb.Hsapiens.UCSC.hg19.knownGene包
NCBI資料庫：ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/

R實現代碼示例

rm(list=ls())na=read.table(choose.files(),sep = t,stringsAsFactors = F,header = T) # 選擇你下的CCDs文件ntmp <- apply(a[1:100,], 1, function(gene){ # 取前100行數據分析調試n# gene=a[3,]n x=gene[10] # Column10 外顯子坐標位置列n if(grepl(],x)){ # 判斷x中是否存在有]這樣的符號，如果有就利用正則替換掉。n x=sub([,,x)n x=sub(],,x)n # 這個時候得到的對象還是像這樣的「880073-880179, 880436-880525……」n tmp <- strsplit(as.character(x),,)[[1]]# 我們先從逗號開始分割成小塊n start <- as.numeric(unlist(lapply(tmp,function(y){# 取開始位點n strsplit(as.character(y),-)[[1]][1]n })))n end <- as.numeric(unlist(lapply(tmp,function(y){ # 取結束位點n strsplit(as.character(y),-)[[1]][2]n })))n gene_d <- data.frame(gene=gene[3], # 將基因名，染色體，開始、結束位點綁定為數據框n chr=gene[1],n start=start,n end=endn )n return (gene_d)#返回數據框n }n}) ntmp_pos=c() # 構造一個空的向量nlapply(tmp[1:10], function(x){ # 取前10個list文件計算調試n# print(x)n if ( !is.null(x)){n apply(x, 1,function(y){n #print(y)n for ( i in as.numeric(y[3]):as.numeric(y[4]) ) # y[3]為坐標起點，y[4]為終止坐標，歷編n tmp_pos <<- c(tmp_pos,paste0(y[2],"-",i))n })nn }n})nlength(tmp_pos) # 計算exon的長度nlength(unique(tmp_pos)) # 計算去重後的exon的長度n

hg19基因組序列的一些探究^3

題目

求：hg19 每條染色體長度，每條染色體N的含量，GC含量。

（高級作業：蛋白編碼區域的GC含量會比基因組其它區域的高嗎？）

測試數據

hg19基因組序列下載

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz # 也可以在瀏覽器上下載ntar xvzf chromFa.tar.gzncat *.fa > hg19.fanrm chr*.fa # 先把著急刪，我待會可能要那他測試運行速度n

簡單測試數據>chr_1
ATCGTCGaaAATGAANccNNttGTA
AGGTCTNAAccAAttGggG
>chr_2
ATCGAATGATCGANNNGccTA
AGGTCTNAAAAGG
>chr_3
ATCGTCGANNNGTAATggGA
AGGTCTNAAAAGG
>chr_4

ATCGTCaaaGANNAATGANGgggTA

PERL代碼示例n

單行命令 {-}

perl -alne {if(/^>/){$chr=$_}else{ $A_count{$chr}+=($_=~tr/Aa//); $T_count{$chr}+=($_=~tr/Tt//);$C_count{$chr}+=($_=~tr/Cc//); $G_count{$chr}+=($_=~tr/Gg//); $N_count{$chr}+=($_=~tr/Nn//); }}END{print "$_t$A_count{$_}t$T_count{$_}t$C_count{$_}t$G_count{$_}t$N_count{$_}" foreach sort keys %N_count} test.fa

完整代碼 {-}

while (<>){nchomp;nif(/^>/){n$chr=$_n}nelse{ n$A_count{$chr}+=($_=~tr/Aa//);n$T_count{$chr}+=($_=~tr/Tt//);n$C_count{$chr}+=($_=~tr/Cc//); n$G_count{$chr}+=($_=~tr/Gg//);n$N_count{$chr}+=($_=~tr/Nn//); n}n} nforeach (sort keys %N_count){n$length = $A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_}+$N_count{$_};n$gc = ($G_count{$_}+$C_count{$_})/($A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_});nprint "$_t$A_count{$_}t$T_count{$_}t$C_count{$_}t$G_count{$_}t$N_count{$_}t$lengtht$gcn" n}n

參考結果{-}

結果如下；

>chr_1 13 10 7 10 4n>chr_2 11 6 5 8 4n>chr_3 10 6 3 10 4n>chr_4 9 4 2 7 3n

hg38每條染色體的基因、轉錄本分布^4

題目

對GTF注釋文件進行探究，統計每條染色體基因數、轉錄本數、內含子數、外顯子數。

（高級作業：下載human/rat/mouse/dog/cat/chicken等物種的gtf注釋文件（http://asia.ensembl.org/info/data/ftp/index.html），編寫函數實現對多個GTF文件進行批量統計染色體基因、轉錄本的分布及轉錄本外顯子個數；繼續探索回答以下問題：所有基因平均有多少個轉錄本？所有轉錄本平均有多個exon和intron？如果要比較多個資料庫呢（gencode/UCSC/NCBI？）？如果把基因分成多個類型呢？protein coding gene，pseudogene，lncRNA還有miRNA的基因？它們的特徵又有哪些變化呢？）

測試數據

wget -c ftp://ftp.ensembl.org/pub/release-87/gtf/homo_sapiens/Homo_sapiens.GRCh38.87.chr.gtf.gzngzip -d Homo_sapiens.GRCh38.87.chr.gtf.gzn

代碼示例

# 每條染色體的基因個數nzcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne {print if $F[2] eq "gene" } |cut -f 1 |sort |uniq -cn# 基因分類nzcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne {next unless $F[2] eq "gene" ;/gene_biotype "(.*?)";/;print $1} |sort |uniq -cn

多個同樣行列式文件的合併^5

題目

將htseq-count生成的所有獨立樣本文件進行合併（每個樣品對應一個文件，包括了所有基因表達量）。

希望通過編程處理每個文件得到輸出的表達矩陣（行名是基因名，列名是樣品名），如下所示：

測試數據

模擬數據

用perl腳本模仿htseq-count計算每個樣本所有的基因表達量的輸出獨立樣本文件：

## 首先新建文件tmp.sh 輸入這個代碼：nperl -le {print "gene_$_t".int(rand(1000)) foreach 1..99}n## 然後用perl腳本調用這個tmp.sh文件：nperl -e system(" bash tmp.sh >$_.txt") foreach a..zn## 這樣就生成了a~z這26個樣本的counts文件n

第一列是基因，第二列是該基因的counts值，共有a~z這26個樣本的counts文件，需要合併成一個大的行列式，這樣才能導入到R裡面做差異分析，如果手工用excel表格做，當然是可以的，但是太麻煩，如果有500個樣本，正常人都不會去手工做了，需要編程。

每個樣本的基因順序並不一致，這時候你應該怎麼做呢？

真實數據

實際需求如下：GSE48213裡面有56個文件，需要合併成一個表達矩陣，來根據cell-line的不同，分組做差異分析。可以查看paper

wget -c ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE48nnn/GSE48213/suppl/GSE48213_RAW.tarntar -xf GSE48213_RAW.tarngzip -d *.gzn

代碼示例

## 首先在GSE48213_RAW目錄裡面生成tmp.txt文件，使用shell腳本：nawk {print FILENAME"t"$0} * |grep -v EnsEMBL_Gene_ID >tmp.txtn## 然後把tmp.txt導入R語言裡面用reshape2處理即可！nsetwd(tmp/GSE48213_RAW/)na=read.table(tmp.txt,sep = t,stringsAsFactors = F)nlibrary(reshape2)nfpkm <- dcast(a,formula = V2~V1)n

根據GTF畫基因的多個轉錄本結構^6

題目

從NCBI,ENSEMBL,UCSC,GENCODE資料庫下載各種GTF注釋文件，編寫代碼得到所有基因的轉錄本個數，以及每個轉錄本的外顯子的坐標，繪製如下轉錄本結構圖：

比如對這個ANXA1基因來說，非常多的轉錄本，但是基因的起始終止坐標，是所有轉錄本起始終止坐標的極大值和極小值。同時，它是一個閉合基因，因為它存在一個轉錄本的起始終止坐標等於該基因的起始終止坐標。可以看到它的外顯子並不是非常整齊的，雖然多個轉錄本會共享某些外顯子，但是也存在某些外顯子比同區域其它外顯子長的現象。

測試數據

wget -c http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/gencode.v7.annotation_goodContig.gtf.gzngzip -d gencode.v7.annotation_goodContig.gtf.gzn

R實現代碼示例

rm(list=ls())nn## http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/gencode.v7.annotation_goodContig.gtf.gznnsetwd(tmp)ngtf <- read.table(gencode.v7.annotation_goodContig.gtf.gz,stringsAsFactors = F,n header = F,comment.char = "#",sep = tn )ntable(gtf[,2])ngtf <- gtf[gtf[,2] ==HAVANA,]ngtf <- gtf[grepl(protein_coding,gtf[,9]),]nnlapply(gtf[1:10,9], function(x){n y=strsplit(x,;) n})nngtf$gene <- lapply(gtf[,9], function(x){n y <- strsplit(x,;)[[1]][5]n strsplit(y,s)[[1]][3]n }n)ndraw_gene = TP53nstructure = gtf[gtf$gene==draw_gene,]nncolnames(structure) =c(n chr,db,record,start,end,tmp1,tmp2,tmp3,tmp4,genen)ngene_start <- min(c(structure$start,structure$end))ngene_end <- max(c(structure$start,structure$end))ntmp_min=min(c(structure$start,structure$end))nstructure$new_start=structure$start-tmp_minnstructure$new_end=structure$end-tmp_minntmp_max=max(c(structure$new_start,structure$new_end))nnum_transcripts=nrow(structure[structure$record==transcript,])ntmp_color=rainbow(num_transcripts)nnx=1:tmp_max;y=rep(num_transcripts,length(x))n#x=10000:17000;y=rep(num_transcripts,length(x))nplot(x,y,type = n,xlab=,ylab = ,ylim = c(0,num_transcripts+1))ntitle(main = draw_gene,sub = paste("chr",structure$chr,":",gene_start,"-",gene_end,sep=""))nj=0;ntmp_legend=c()nfor (i in 1:nrow(structure)){n tmp=structure[i,]n if(tmp$record == transcript){n j=j+1n tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))n }n if(tmp$record == exon) lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)n}n

參考結果

下載最新版的KEGG信息，並且解析好^7

題目

下載最新版的KEGG注釋文本文件，編寫腳本整理成kegg的pathway的ID與基因ID的對應格式。

測試數據

1 首先打開KEGG官方網站，網頁中展示出了各個物種的分類、拉丁名稱、英文名稱等信息。

2 直接網頁中搜索（Ctrl + F）需要下載的物種英文名稱或拉丁名。如果不確定物種名稱，網站中提供了詳細的分類系統，也可根據前面的物種分類信息進行查找。
本文以擬南芥為例，搜索「Arabidopsis thaliana」即可找到。找到後點擊物種名稱前的3個字母縮寫鏈接（下圖紅色框中的位置）。

3 進入後的網頁中包含了物種的一些基因組信息，點擊上方的「Brite hierarchy」，進入後再點擊「KEGG Orthology (KO)」；

4 在跳轉出的網頁中點擊「Download htext」，彈出下載窗口進行下載，國外網站有時會出現無法下載的情況，多試幾次即可；

5 當然，下載好之後還沒有結束。下載得到文本文件，可以看到裡面的結構層次非常清楚，C開頭的就是kegg的pathway的ID所在行，D開頭的就是屬於它的kegg的所有的基因。A,B是kegg的分類，總共是6個大類，42個小類。

PERL代碼示例

perl -alne {if(/^C/){/PATH:hsa(d+)/;$kegg=$1}else{print "$keggt$F[1]" if /^D/ and $kegg;}} hsa00001.keg >kegg2gene.txtn

參考結果

寫超幾何分布檢驗^8

題目

學習GO/KEGG的富集分析的原理，編寫代碼實現超幾何分布檢驗，將得到的結果與測試數據中的kegg.enrichment.html進行比較。

(P值的計算：C(k,M)*C(n-k,N-M)/C(n,N) )

測試數據

kegg2gene（第六講kegg數據解析結果）

暫時不用最新版的kegg注釋數據，為了能夠統一答案

差異基因list和背景基因list(R代碼)

source("http://bioconductor.org/biocLite.R")nbiocLite("org.Hs.eg.db")nbiocLite("KEGG.db")nbiocLite("GOstats")nbiocLite("hgu95av2.db")nnlibrary(org.Hs.eg.db)nlibrary(KEGG.db)nlibrary(GOstats)nlibrary("hgu95av2.db")nn##得到kegg2gene.list(KEGG注釋信息)nntmp=toTable(org.Hs.egPATH)nwrite.table(tmp,kegg2gene.list.txt,quote = F,row.names = F)n## 得到universeGeneIds.txt(背景基因list)nls(package:hgu95av2.db)nuniverseGeneIds <- unique(mappedRkeys(hgu95av2ENTREZID))nwrite.table(universeGeneIds,universeGeneIds.txt,quote = F,row.names = F)nn## 得到diff_gene.txt(差異基因list)nnset.seed(123456.789)ndiff_gene = sample(universeGeneIds,300)nwrite.table(diff_gene,diff_gene.txt,quote = F,row.names = F)nn## 得到kegg.enrichment.html(富集分析結果)nnannotationPKG=org.Hs.eg.dbnhyperG.params = new("KEGGHyperGParams", geneIds=diff_gene, universeGeneIds=universeGeneIds, annotation=annotationPKG,ncategoryName="KEGG", pvalueCutoff=1, testDirection = "over")nKEGG.hyperG.results = hyperGTest(hyperG.params);nhtmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))n

R代碼示例

下載鏈接： https://github.com/jmzeng1314/humanid/blob/master/R/hyperGtest_jimmy.R

library("hgu95av2.db")[/align][align=left]ls(package:hgu95av2.db)nuniverseGeneIds <- unique(mappedRkeys(hgu95av2ENTREZID))nset.seed(123456.789)ndiff_gene = sample(universeGeneIds,300)nlibrary(org.Hs.eg.db)nlibrary(KEGG.db)ntmp=toTable(org.Hs.egPATH)nGeneID2kegg<<- tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)nkegg2GeneID<<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)n#results <- hyperGtest_jimmy(GeneID2kegg,kegg2GeneID,diff_gene,universeGeneIds)nnGeneID2Path=GeneID2keggnPath2GeneID=kegg2GeneIDnndiff_gene_has_path=intersect(diff_gene,names(GeneID2Path))nn=length(diff_gene) #306nN=length(universeGeneIds) #5870nresults=c()nnfor (i in names(Path2GeneID)){n #i="04672"n M=length(intersect(Path2GeneID[[i]],universeGeneIds))n #print(M)nn if(M<5)n nextn exp_count=n*M/Nn #print(paste(n,N,M,sep="t"))n k=0n for (j in diff_gene_has_path){n if (i %in% GeneID2Path[[j]]) k=k+1n }n OddsRatio=k/exp_countn if (k==0) nextn p=phyper(k-1,M, N-M, n, lower.tail=F)n #print(paste(i,p,OddsRatio,exp_count,k,M,sep=" "))n results=rbind(results,c(i,p,OddsRatio,exp_count,k,M))n}ncolnames(results)=c("PathwayID","Pvalue","OddsRatio","ExpCount","Count","Size")nresults=as.data.frame(results,stringsAsFactors = F)nresults$p.adjust = p.adjust(results$Pvalue,method = BH)nresults=results[order(results$Pvalue),]nrownames(results)=1:nrow(results)n

ID轉換^9

題目

probe_id，gene_id，gene_name, symbol之間的轉換。

測試數據

需要轉換的探針ID(變數my_probe)

rm(list=ls())nlibrary("hgu95av2.db")nls(package:hgu95av2.db)nprobe2entrezID=toTable(hgu95av2ENTREZID)nprobe2symbol=toTable(hgu95av2SYMBOL)nprobe2genename=toTable(hgu95av2GENENAME)nnmy_probe = sample(unique(mappedLkeys(hgu95av2ENTREZID)),30)nntmp1 = probe2symbol[match(my_probe,probe2symbol$probe_id),]ntmp2 = probe2entrezID[match(my_probe,probe2entrezID$probe_id),]ntmp3 = probe2genename[match(my_probe,probe2genename$probe_id),]nnwrite.table(my_probe,my_probe.txt,quote = F,col.names = F,row.names =F)nwrite.table(tmp1$symbol,my_symbol.txt,quote = F,col.names = F,row.names =F)nwrite.table(tmp2$gene_id,my_geneID.txt,quote = F,col.names = F,row.names =F)n

下載對應關係
ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
illumina晶元的探針

所有bioconductor支持的晶元平台對應關係：通過bioconductor包來獲取所有的晶元探針與gene的對應關係；可以從NCBI的GPL平台下載：http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL6947；也可以直接載入對應的包；或者直接去公司的主頁下載manifest文件。

library("illuminaHumanv4.db")nls(package:illuminaHumanv4.db)nnprobe2entrezID=toTable(illuminaHumanv4ENTREZID)nprobe2symbol=toTable(illuminaHumanv4SYMBOL)nprobe2genename=toTable(illuminaHumanv4GENENAME)nnmy_probe = sample(unique(mappedLkeys(illuminaHumanv4ENTREZID)),30)nnprobe2symbol[match(my_probe,probe2symbol$probe_id),]nprobe2entrezID[match(my_probe,probe2entrezID$probe_id),]nprobe2genename[match(my_probe,probe2genename$probe_id),]n

R代碼示例

基因的轉換：運行下面的R代碼，得到的my_symbol_gene和my_entrez_gene就是需要轉換的ID。

library("illuminaHumanv4.db")nls(package:illuminaHumanv4.db)nmy_entrez_gene = sample(unique(mappedRkeys(illuminaHumanv4ENTREZID)),30)nmy_symbol_gene = sample(unique(mappedRkeys(illuminaHumanv4SYMBOL)),30)nnlibrary("org.Hs.eg.db")nls(package:org.Hs.eg.db)nnentrezID2symbol <- toTable(org.Hs.egSYMBOL)nnentrezID2symbol[match(my_entrez_gene,entrezID2symbol$gene_id),]nentrezID2symbol[match(my_symbol_gene,entrezID2symbol$symbol),]n

根據指定染色體及坐標得到序列^10

題目

參考基因組hg19，指定染色體及坐標」chr5」,」8397384」，編寫程序得到這個坐標以及它上下一個鹼基。（機器無法計算hg19，則使用測試數據，指定坐標是 3號染色體的第6個鹼基。）

測試數據

>chr_1nATCGTCGaaAATGAANccNNttGTAnAGGTCTNAAccAAttGggGn>chr_2nATCGAATGATCGANNNGccTAnAGGTCTNAAAAGGn>chr_3nATCGTCGANNNGTAATggGAnAGGTCTNAAAAGGn>chr_4nATCGTCaaaGANNAATGANGgggTAn

根據指定染色體及坐標得到位置信息^11

題目

基因的chr,start,end都是已知的（坐標是hg38系統），任意給定基因組的chr:pos(chr1:2075000-2930999), 判斷該區間在哪個基因上面？（可用現成軟體bedtools）

測試數據

chr7 148697841 148698941nchr7 148698942 148699029nchr7 148699911 148701053nchr7 148701109 148701307nchr7 148701354 148702694nchr7 148703100 148703520nchr7 148703831 148704175nchr7 148704484 148704734nchr7 148704857 148705937nchr7 148706271 148706671n

把文件內容按照染色體分開寫出^12

題目

根據染色體把文件拆分成1~22和其它染色體的兩個文件呢？

測試數據

wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txtn

把文件改成bed格式，如下：

chr2 43995310 43995986nchr17 49788603 49789067nchr17 59565573 59566163nchr19 8390308 8390745nchr12 49188033 49189033nchr7 974903 975570nchr7 98878532 98879500nchr7 44044672 44045322nchr1 153634052 153634772nchr11 60905850 60906575n

PERL代碼示例

use FileHandle;nforeach( 1..22 ){n $normal_chr{"chr".$_}=1 ;n open $fh{"chr".$_},">chr$_.txt" or die;n}nopen other,">chr_other.txt" or die;nopen FH,a.bed;nwhile(<FH>){n chomp;n @F=split;n if(exists $normal_chr{$F[0]}){n $fh{$F[0]}->print( "$_n" );n }else{ n print other $_;n }n}nforeach( 1..22 ){ n close $fh{$_};n}nclose other;n

JSON格式數據的格式化^13

題目

學習json格式，下載測試數據，從該json文件裡面提取：technique factor target principal_investigator submission label category type Developmental-Stage organism key這幾列信息。

測試數據

http://biotrainee.com/jbrowse/JBrowse-1.12.1/sample_data/json/modencode/modencodeMetaData.json

PERL代碼示例

#!/usr/bin/env perlnuse strict;nuse warnings;nuse autodie :all;nuse 5.10.0;nnuse JSON 2;nnmy $data = from_json( do { local $/; open my $f, <, $ARGV[0]; scalar <$f> } );nnmy @fields = qw( technique factor target principal_investigator submission label category type Developmental-Stage organism key );nnsay join ,, map ""$_"", @fields;nnfor my $item ( @{$data->{items}} ) {n $item->{key} = $item->{label};n no warnings uninitialized;n for my $track ( @{$item->{Tracks}} ) {n $item->{label} = $track;n say join ,, map ""$_"", @{$item}{@fields};n }n}n

參考結果

完成之後的結果應該是：http://biotrainee.com/jbrowse/JBrowse-1.12.1/sample_data/json/modencode/modencodeMetaData.csv

多個探針對應一個基因，取平均值、最大值 [^14]

[^14]: 多個探針對應一個基因，取平均值、最大值

題目

編寫腳本對多個探針對應一個基因，取平均值、最大值。

測試數據

R代碼示例

# 平均值nBiocInstaller::biocLite(CLL)nBiocInstaller::biocLite(hgu95av2.db)nnlibrary(hgu95av2.db)nlibrary(CLL)ndata(sCLLex)nsCLLex=sCLLex[,1:8] ## 樣本太多，我就取前面8個ngroup_list=sCLLex$DiseasenexprSet=exprs(sCLLex)nexprSet=as.data.frame(exprSet)nexprSet$probe_id=rownames(exprSet)nhead(exprSet)nprobe2symbol=toTable(hgu95av2SYMBOL)ndat=merge(exprSet,probe2symbol,by=probe_id)nresults=t(sapply(split(dat,dat$symbol),function(x) colMeans(x[,1:(ncol(x)-1)])))n

把counts矩陣轉換成RPKM矩陣[^15]

[^15]: 把counts矩陣轉換成RPKM矩陣

題目

編寫腳本將hisat2+htseq-counts之後得到reads的counts矩陣轉換成RPKM矩陣。

測試數據

wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_mouse/CCDS.current.txtngrep -v ^# CCDS.current.txt | perl -alne {/[(.*?)]/;$len=0;foreach(split/,/,$1){@tmp=split/-/;$len+=($tmp[1]-$tmp[0])};$h{$F[2]}=$len if $len >$h{$F[2]}} END{print "$_t$h{$_}" foreach sort keys %h} >mm10_ccds_length.txtn

R代碼示例

genes_len=read.table("mm10_ccds_length.txt",stringsAsFactors=F)nhead(genes_len)n V1 V2n1 -343C11.2 1139n2 00R_AC107638.2 1060n3 00R_Pgap2 1676n4 0610005C13Rik 7363n5 0610006L08Rik 34995n6 0610007P14Rik 9074nncolnames(genes_len)<- c("GeneName","Len")nnhead(exprSet)n GSM860181_priSG-A GSM860182_SG-A GSM860183_SG-B GSM860184_lepSCn00R_AC107638.2 0 1 0 1n0610005C13Rik 20 22 11 27n0610006L08Rik 0 0 0 2n0610007P14Rik 2075 1785 1269 1430n0610009B22Rik 256 376 300 386n0610009E02Rik 22 22 16 28nnexprSet<-exprSet[ rownames(exprSet) %in% genes_len$GeneName ,]nntotal_count<- colSums(exprSet)nnneededGeneLength=genes_len[ match(rownames(exprSet), genes_len$GeneName) ,2] nnrpkm <- t(do.call( rbind,lapply(1:length(total_count),function(i){n 10^9*exprSet[,i]/neededGeneLength/total_count[i]n}) ))nhead(rpkm)nrownames(rpkm)=rownames(exprSet)ncolnames(rpkm)=colnames(exprSet)nnwrite.table(rpkm,file="rpkm.txt",sep="t",quote=F)nnnlibrary(TxDb.Mmusculus.UCSC.mm10.knownGene)ntxdb=TxDb.Mmusculus.UCSC.mm10.knownGenengt=transcriptsBy(txdb,by="gene")nlapply(gt[1:40],function(x) max(width(x)))nlibrary(org.Mm.eg.db)nnmykeys=ncolumns(txdb);keytypes(txdb)nneededcols <- c("GENEID", "TXSTRAND", "TXCHROM")nselect(txdb, keys=mykeys, columns=neededcols, keytype="TXNAME")n

參考結果

對有臨床信息的表達矩陣批量做生存分析^16

題目

使用R實現生存分析：

用my.surv <- surv(OS_MONTHS,OS_STATUS==DECEASED)構建生存曲線；

用kmfit2 <- survfit(my.surv~TUMOR_STAGE_2009)來做某一個因子的KM生存曲線；

用survdiff(my.surv~type, data=dat)來看看這個因子的不同水平是否有顯著差異，其中默認用是的logrank test 方法；

用coxph(Surv(time, status) ~ ph.ecog + tt(age), data=lung)來檢測自己感興趣的因子是否受其它因子(age,gender等等)的影響。

代碼示例

rm(list=ls())n## 50 patients and 200 genes ndat=matrix(rnorm(10000,mean=8,sd=4),nrow = 200)nrownames(dat)=paste0(gene_,1:nrow(dat))ncolnames(dat)=paste0(sample_,1:ncol(dat))nos_years=abs(floor(rnorm(ncol(dat),mean = 50,sd=20)))nos_status=sample(rep(c(LIVING,DECEASED),25))nnlibrary(survival)nmy.surv <- Surv( os_years,os_status==DECEASED)n## The status indicator, normally 0=alive, 1=dead. Other choices are TRUE/FALSE (TRUE = death) or 1/2 (2=death). n## And most of the time we just care about the time od DECEASED .nnfit.KM=survfit(my.surv~1)nfit.KMnplot(fit.KM)nnlog_rank_p <- apply(dat, 1, function(values1){n group=ifelse(values1>median(values1),high,low)n kmfit2 <- survfit(my.surv~group)n #plot(kmfit2)n data.survdiff=survdiff(my.surv~group)n p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)n})nnames(log_rank_p[log_rank_p<0.05])nngender= ifelse(rnorm(ncol(dat))>1,male,female)nage=abs(floor(rnorm(ncol(dat),mean = 50,sd=20)))n## gender and age are similar with group(by gene expression)nncox_results <- apply(dat , 1, function(values1){n group=ifelse(values1>median(values1),high,low)n survival_dat <- data.frame(group=group,gender=gender,age=age,stringsAsFactors = F)n m=coxph(my.surv ~ age + gender + group, data = survival_dat)nn beta <- coef(m)n se <- sqrt(diag(vcov(m)))n HR <- exp(beta)n HRse <- HR * senn #summary(m)n tmp <- round(cbind(coef = beta, se = se, z = beta/se, p = 1 - pchisq((beta/se)^2, 1),n HR = HR, HRse = HRse,n HRz = (HR - 1) / HRse, HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),n HRCILL = exp(beta - qnorm(.975, 0, 1) * se),n HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)n return(tmp[grouplow,])nn})ncox_results[,cox_results[4,]<0.05]n

PS: 裡面的os_years應該修改為os_month，因為總的生存期不應該長達幾十年，而且因為ag和os_years都是隨機生成的，可能會出現很不符合自然科學的現象。但是這個是模擬數據，請大家不用較真。

對多個差異分析結果直接取交集並集[^17]

[^17]: 對多個差異分析結果直接取交集並集

題目

編寫腳本對每兩個差異分析結果計算基因交集個數與基因並集個數的比值，得到一個比值矩陣。

測試數據

用perl單行命令模擬數據：

perl -e BEGIN{use List::Util qw/shuffle/;@gene=A..Z}{foreach(1..10){@this_genes=@gene[(shuffle(0..25))[0..int(rand(20))+4]];print "comparison_$_ <- ";print join(";",@this_genes);print "n"}}n

總共10次差異分析，每次差異分析得到的差異基因在同一行的後面用大小字母表示。

comparison_1 -> I;G;E;V;C;K;B;Wncomparison_2 -> G;E;N;H;Y;M;L;S;K;A;J;O;D;P;R;U;Q;F;Z;Cncomparison_3 -> Y;V;U;N;H;K;I;P;S;F;D;X;G;C;Z;J;Q;T;W;O;E;Mncomparison_4 -> N;T;K;B;H;Z;W;C;Q;I;V;F;D;S;R;Y;J;X;P;O;G;L;Ancomparison_5 -> G;J;A;H;W;T;Z;E;Y;S;Rncomparison_6 -> Z;M;D;R;P;G;L;W;Y;U;X;E;A;S;T;I;Hncomparison_7 -> H;Z;T;O;W;Q;M;X;J;N;U;K;F;P;I;C;S;Y;A;Bncomparison_8 -> A;R;L;T;W;Q;S;F;P;X;E;V;Y;G;K;J;Z;Cncomparison_9 -> J;X;K;D;O;H;L;F;C;P;R;Nncomparison_10 -> G;S;K;H;C;O;W;F;Q;Xn

R代碼示例

a=readLines(test.txt)nn=unlist(lapply(a , function(x){n trimws(strsplit(x,<-)[[1]][1])n}))nv=lapply(a , function(x){n trimws(strsplit(trimws(strsplit(x,<-)[[1]][2]),;)[[1]])n})nnames(v)=nntmp=unlist(lapply(v, function(x){n lapply(v, function(y){n p1=length(intersect(x,y))n p2=length(union(x,y))n return(p1/p2)n })n}))nout=matrix(tmp,nrow = length(n))nrownames(out)=nncolnames(out)=nnout[1:5,1:5]n

參考結果

結果的前五行如下：

comparison_1 comparison_2 comparison_3 comparison_4 comparison_5ncomparison_1 1.0000000 0.1666667 0.3043478 0.2916667 0.1875000ncomparison_2 0.1666667 1.0000000 0.6800000 0.6538462 0.4090909ncomparison_3 0.3043478 0.6800000 1.0000000 0.7307692 0.3750000ncomparison_4 0.2916667 0.6538462 0.7307692 1.0000000 0.4166667ncomparison_5 0.1875000 0.4090909 0.3750000 0.4166667 1.0000000n

根據GTF格式的基因注釋文件得到人所有基因的染色體坐標[^18]

[^18]: 根據GTF格式的基因注釋文件得到人所有基因的染色體坐標

題目

從gencode資料庫裡面可以下載所有的gtf文件，編寫腳本得到基因的染色體、起始終止坐標如下：

[jianmingzeng@gencode]$ head protein_coding.hg19.position nchr1 69091 70008 OR4F5nchr1 367640 368634 OR4F29nchr1 621096 622034 OR4F16nchr1 859308 879961 SAMD11nchr1 879584 894689 NOC2Lnchr1 895967 901095 KLHL17nchr1 901877 911245 PLEKHN1nchr1 910584 917473 PERM1nchr1 934342 935552 HES4nchr1 936518 949921 ISG15n[jianmingzeng@gencode]$ head protein_coding.hg38.position nchr1 69091 70008 OR4F5nchr1 182393 184158 FO538757.2nchr1 184923 200322 FO538757.1nchr1 450740 451678 OR4F29nchr1 685716 686654 OR4F16nchr1 923928 944581 SAMD11nchr1 944204 959309 NOC2Lnchr1 960587 965715 KLHL17nchr1 966497 975865 PLEKHN1nchr1 975204 982093 PERM1n