生物資訊程式設計實戰題
阿新 • • 發佈：2018-11-10
目錄
18.根據GTF格式的基因註釋檔案得到人所有基因的染色體座標
1.生信程式設計很簡單
程式語言系統入門
題目
對FASTQ的操作
- 5,3段截掉幾個鹼基
- 序列長度分佈統計
- FASTQ 轉換成 FASTA
- 統計鹼基個數及GC%
對FASTA的操作
- 取互補序列
- 取反向序列
- DNA to RNA
- 大小寫字母形式輸出
- 每行指定長度輸出序列
- 按照序列長度/名字排序
- 提取指定ID的序列
- 隨機抽取序列
高階難度
- 根據座標取序列
- 多檔案合併
- 根據ID列表取序列
- GTF檔案探索
- 簡併鹼基的引物序列還原成多條序列
- snp進行註釋並格式化輸出
下載安裝bowtie2（內含測試資料）
cd ~/biosoft
mkdir bowtie && cd bowtie
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip
unzip bowtie2-2.2.9-linux-x86_64.zip
2.人類基因組的外顯子區域的長度
題目
下載人類外顯子的座標檔案，編寫程式碼統計外顯子區域的長度。
測試資料
- Rbioconductor的TxDb.Hsapiens.UCSC.hg19.knownGene包
- NCBI資料庫：ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/
R實現程式碼示例
rm(list=ls())
a=read.table(choose.files(),sep = '\t',stringsAsFactors = F,header = T) # 選擇你下的CCDs檔案
tmp <- apply(a[1:100,], 1, function(gene){ # 取前100行資料分析除錯
# gene=a[3,]
x=gene[10] # Column10 外顯子座標位置列
if(grepl('\\]',x)){ # 判斷x中是否存在有]這樣的符號，如果有就利用正則替換掉。
x=sub('\\[','',x)
x=sub('\\]','',x)
# 這個時候得到的物件還是像這樣的“880073-880179, 880436-880525……”
tmp <- strsplit(as.character(x),',')[[1]]# 我們先從逗號開始分割成小塊
start <- as.numeric(unlist(lapply(tmp,function(y){# 取開始位點
strsplit(as.character(y),'-')[[1]][1]
})))
end <- as.numeric(unlist(lapply(tmp,function(y){ # 取結束位點
strsplit(as.character(y),'-')[[1]][2]
})))
gene_d <- data.frame(gene=gene[3], # 將基因名，染色體，開始、結束位點繫結為資料框
chr=gene[1],
start=start,
end=end
)
return (gene_d)#返回資料框
}
})
tmp_pos=c() # 構造一個空的向量
lapply(tmp[1:10], function(x){ # 取前10個list檔案計算除錯
# print(x)
if ( !is.null(x)){
apply(x, 1,function(y){
#print(y)
for ( i in as.numeric(y[3]):as.numeric(y[4]) ) # y[3]為座標起點，y[4]為終止座標，歷編
tmp_pos <<- c(tmp_pos,paste0(y[2],"-",i))
})
}
})
length(tmp_pos) # 計算exon的長度
length(unique(tmp_pos)) # 計算去重後的exon的長度
3.hg19基因組序列的一些探究
題目
求：hg19 每條染色體長度，每條染色體N的含量，GC含量。 （高階作業：蛋白編碼區域的GC含量會比基因組其它區域的高嗎？ ）
測試資料
- hg19基因組序列下載
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz # 也可以在瀏覽器上下載
tar xvzf chromFa.tar.gz
cat *.fa > hg19.fa
rm chr*.fa # 先把著急刪，我待會可能要那他測試執行速度
- 簡單測試資料
>chr_1
ATCGTCGaaAATGAANccNNttGTA
AGGTCTNAAccAAttGggG
>chr_2
ATCGAATGATCGANNNGccTA
AGGTCTNAAAAGG
>chr_3
ATCGTCGANNNGTAATggGA
AGGTCTNAAAAGG
>chr_4
ATCGTCaaaGANNAATGANGgggTA
Perl程式碼示例
單行命令
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 (<>){
chomp;
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//);
}
}
foreach (sort keys %N_count){
$length = $A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_}+$N_count{$_};
$gc = ($G_count{$_}+$C_count{$_})/($A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_});
print "$_\t$A_count{$_}\t$T_count{$_}\t$C_count{$_}\t$G_count{$_}\t$N_count{$_}\t$length\t$gc\n"
}
參考結果{-}
結果如下；
>chr_1 13 10 7 10 4
>chr_2 11 6 5 8 4
>chr_3 10 6 3 10 4
>chr_4 9 4 2 7 3
4.hg38每條染色體的基因、轉錄本分佈
題目
對GTF註釋檔案進行探究，統計每條染色體基因數、轉錄本數、內含子數、外顯子數。
高階作業：去ENSEMBL資料庫 下載human/rat/mouse/dog/cat/chicken等物種的gtf註釋檔案編寫函式實現對多個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.gz
gzip -d Homo_sapiens.GRCh38.87.chr.gtf.gz
程式碼示例
# 每條染色體的基因個數
zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{print if $F[2] eq "gene" }' |cut -f 1 |sort |uniq -c
# 基因分類
zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{next unless $F[2] eq "gene" ;/gene_biotype "(.*?)";/;print $1}' |sort |uniq -c
5.多個同樣行列式檔案的合併
題目
將htseq-count生成的所有獨立樣本檔案進行合併（每個樣品對應一個檔案，包括了所有基因表達量）。 希望通過程式設計處理每個檔案得到輸出的表達矩陣（行名是基因名，列名是樣品名），如下所示：
模擬資料
用perl指令碼模仿htseq-count計算每個樣本所有的基因表達量的輸出獨立樣本檔案：
## 首先新建檔案tmp.sh 輸入這個程式碼：
perl -le '{print "gene_$_\t".int(rand(1000)) foreach 1..99}'
## 然後用perl指令碼呼叫這個tmp.sh檔案：
perl -e 'system(" bash tmp.sh >$_.txt") foreach a..z'
## 這樣就生成了a~z這26個樣本的counts檔案
第一列是基因，第二列是該基因的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.tar
tar -xf GSE48213_RAW.tar
gzip -d *.gz
程式碼示例
## 首先在GSE48213_RAW目錄裡面生成tmp.txt檔案，使用shell指令碼：
awk '{print FILENAME"\t"$0}' * |grep -v EnsEMBL_Gene_ID >tmp.txt
## 然後把tmp.txt匯入R語言裡面用reshape2處理即可！
setwd('tmp/GSE48213_RAW/')
a=read.table('tmp.txt',sep = '\t',stringsAsFactors = F)
library(reshape2)
fpkm <- dcast(a,formula = V2~V1)
6.根據GTF畫基因的多個轉錄本結構
題目
從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.gz
gzip -d gencode.v7.annotation_goodContig.gtf.gz
R實現程式碼示例
rm(list=ls())
## http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/gencode.v7.annotation_goodContig.gtf.gz
setwd('tmp')
gtf <- read.table('gencode.v7.annotation_goodContig.gtf.gz',stringsAsFactors = F,
header = F,comment.char = "#",sep = '\t'
)
table(gtf[,2])
gtf <- gtf[gtf[,2] =='HAVANA',]
gtf <- gtf[grepl('protein_coding',gtf[,9]),]
lapply(gtf[1:10,9], function(x){
y=strsplit(x,';')
})
gtf$gene <- lapply(gtf[,9], function(x){
y <- strsplit(x,';')[[1]][5]
strsplit(y,'\\s')[[1]][3]
}
)
draw_gene = 'TP53'
structure = gtf[gtf$gene==draw_gene,]
colnames(structure) =c(
'chr','db','record','start','end','tmp1','tmp2','tmp3','tmp4','gene'
)
gene_start <- min(c(structure$start,structure$end))
gene_end <- max(c(structure$start,structure$end))
tmp_min=min(c(structure$start,structure$end))
structure$new_start=structure$start-tmp_min
structure$new_end=structure$end-tmp_min
tmp_max=max(c(structure$new_start,structure$new_end))
num_transcripts=nrow(structure[structure$record=='transcript',])
tmp_color=rainbow(num_transcripts)
x=1:tmp_max;y=rep(num_transcripts,length(x))
#x=10000:17000;y=rep(num_transcripts,length(x))
plot(x,y,type = 'n',xlab='',ylab = '',ylim = c(0,num_transcripts+1))
title(main = draw_gene,sub = paste("chr",structure$chr,":",gene_start,"-",gene_end,sep=""))
j=0;
tmp_legend=c()
for (i in 1:nrow(structure)){
tmp=structure[i,]
if(tmp$record == 'transcript'){
j=j+1
tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))
}
if(tmp$record == 'exon') lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)
}
參考結果
7.下載最新版的KEGG資訊，並且解析好
題目
下載最新版的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 "$kegg\t$F[1]" if /^D/ and $kegg;}}' hsa00001.keg >kegg2gene.txt
參考結果
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")
biocLite("org.Hs.eg.db")
biocLite("KEGG.db")
biocLite("GOstats")
biocLite("hgu95av2.db")
library(org.Hs.eg.db)
library(KEGG.db)
library(GOstats)
library("hgu95av2.db")
##得到kegg2gene.list(KEGG註釋資訊)
tmp=toTable(org.Hs.egPATH)
write.table(tmp,'kegg2gene.list.txt',quote = F,row.names = F)
## 得到universeGeneIds.txt(背景基因list)
ls('package:hgu95av2.db')
universeGeneIds <- unique(mappedRkeys(hgu95av2ENTREZID))
write.table(universeGeneIds,'universeGeneIds.txt',quote = F,row.names = F)
## 得到diff_gene.txt(差異基因list)
set.seed('123456.789')
diff_gene = sample(universeGeneIds,300)
write.table(diff_gene,'diff_gene.txt',quote = F,row.names = F)
## 得到kegg.enrichment.html(富集分析結果)
annotationPKG='org.Hs.eg.db'
hyperG.params = new("KEGGHyperGParams", geneIds=diff_gene, universeGeneIds=universeGeneIds, annotation=annotationPKG,
categoryName="KEGG", pvalueCutoff=1, testDirection = "over")
KEGG.hyperG.results = hyperGTest(hyperG.params);
htmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))
R程式碼示例
下載連結： https://github.com/jmzeng1314/humanid/blob/master/R/hyperGtest_jimmy.R
library("hgu95av2.db")[/align][align=left]ls('package:hgu95av2.db')
universeGeneIds <- unique(mappedRkeys(hgu95av2ENTREZID))
set.seed('123456.789')
diff_gene = sample(universeGeneIds,300)
library(org.Hs.eg.db)
library(KEGG.db)
tmp=toTable(org.Hs.egPATH)
GeneID2kegg<<- tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)
kegg2GeneID<<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
#results <- hyperGtest_jimmy(GeneID2kegg,kegg2GeneID,diff_gene,universeGeneIds)
GeneID2Path=GeneID2kegg
Path2GeneID=kegg2GeneID
diff_gene_has_path=intersect(diff_gene,names(GeneID2Path))
n=length(diff_gene) #306
N=length(universeGeneIds) #5870
results=c()
for (i in names(Path2GeneID)){
#i="04672"
M=length(intersect(Path2GeneID[[i]],universeGeneIds))
#print(M)
if(M<5)
next
exp_count=n*M/N
#print(paste(n,N,M,sep="\t"))
k=0
for (j in diff_gene_has_path){
if (i %in% GeneID2Path[[j]]) k=k+1
}
OddsRatio=k/exp_count
if (k==0) next
p=phyper(k-1,M, N-M, n, lower.tail=F)
#print(paste(i,p,OddsRatio,exp_count,k,M,sep=" "))
results=rbind(results,c(i,p,OddsRatio,exp_count,k,M))
}
colnames(results)=c("PathwayID","Pvalue","OddsRatio","ExpCount","Count","Size")
results=as.data.frame(results,stringsAsFactors = F)
results$p.adjust = p.adjust(results$Pvalue,method = 'BH')
results=results[order(results$Pvalue),]
rownames(results)=1:nrow(results)
9.ID轉換
題目
probe_id，gene_id，gene_name, symbol之間的轉換。
測試資料
- 需要轉換的探針ID(變數my_probe)
rm(list=ls())
library("hgu95av2.db")
ls('package:hgu95av2.db')
probe2entrezID=toTable(hgu95av2ENTREZID)
probe2symbol=toTable(hgu95av2SYMBOL)
probe2genename=toTable(hgu95av2GENENAME)
my_probe = sample(unique(mappedLkeys(hgu95av2ENTREZID)),30)
tmp1 = probe2symbol[match(my_probe,probe2symbol$probe_id),]
tmp2 = probe2entrezID[match(my_probe,probe2entrezID$probe_id),]
tmp3 = probe2genename[match(my_probe,probe2genename$probe_id),]
write.table(my_probe,'my_probe.txt',quote = F,col.names = F,row.names =F)
write.table(tmp1$symbol,'my_symbol.txt',quote = F,col.names = F,row.names =F)
write.table(tmp2$gene_id,'my_geneID.txt',quote = F,col.names = F,row.names =F)
- 下載對應關係 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")
ls('package:illuminaHumanv4.db')
probe2entrezID=toTable(illuminaHumanv4ENTREZID)
probe2symbol=toTable(illuminaHumanv4SYMBOL)
probe2genename=toTable(illuminaHumanv4GENENAME)
my_probe = sample(unique(mappedLkeys(illuminaHumanv4ENTREZID)),30)
probe2symbol[match(my_probe,probe2symbol$probe_id),]
probe2entrezID[match(my_probe,probe2entrezID$probe_id),]
probe2genename[match(my_probe,probe2genename$probe_id),]
R程式碼示例
基因的轉換：執行下面的R程式碼，得到的my_symbol_gene和my_entrez_gene就是需要轉換的ID。
library("illuminaHumanv4.db")
ls('package:illuminaHumanv4.db')
my_entrez_gene = sample(unique(mappedRkeys(illuminaHumanv4ENTREZID)),30)
my_symbol_gene = sample(unique(mappedRkeys(illuminaHumanv4SYMBOL)),30)
library("org.Hs.eg.db")
ls('package:org.Hs.eg.db')
entrezID2symbol <- toTable(org.Hs.egSYMBOL)
entrezID2symbol[match(my_entrez_gene,entrezID2symbol$gene_id),]
entrezID2symbol[match(my_symbol_gene,entrezID2symbol$symbol),]
10.根據指定染色體及座標得到序列
題目
參考基因組hg19，指定染色體及座標"chr5","8397384"，編寫程式得到這個座標以及它上下一個鹼基。（機器無法計算hg19，則使用測試資料，指定座標是 3號染色體的第6個鹼基。）
測試資料
>chr_1
ATCGTCGaaAATGAANccNNttGTA
AGGTCTNAAccAAttGggG
>chr_2
ATCGAATGATCGANNNGccTA
AGGTCTNAAAAGG
>chr_3
ATCGTCGANNNGTAATggGA
AGGTCTNAAAAGG
>chr_4
ATCGTCaaaGANNAATGANGgggTA
11.根據指定染色體及座標得到位置資訊
題目
基因的chr,start,end都是已知的（座標是hg38系統），任意給定基因組的chr:pos(chr1:2075000-2930999), 判斷該區間在哪個基因上面？（可用現成軟體bedtools）
測試資料
chr7 148697841 148698941
chr7 148698942 148699029
chr7 148699911 148701053
chr7 148701109 148701307
chr7 148701354 148702694
chr7 148703100 148703520
chr7 148703831 148704175
chr7 148704484 148704734
chr7 148704857 148705937
chr7 148706271 148706671
12.把檔案內容按照染色體分開寫出
題目
根據染色體把檔案拆分成1~22和其它染色體的兩個檔案呢？
測試資料
wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/CCDS.current.txt
把檔案改成bed格式，如下：
chr2 43995310 43995986
chr17 49788603 49789067
chr17 59565573 59566163
chr19 8390308 8390745
chr12 49188033 49189033
chr7 974903 975570
chr7 98878532 98879500
chr7 44044672 44045322
chr1 153634052 153634772
chr11 60905850 60906575
Perl程式碼示例
use FileHandle;
foreach( 1..22 ){
$normal_chr{"chr".$_}=1 ;
open $fh{"chr".$_},">chr$_.txt" or die;
}
open other,">chr_other.txt" or die;
open FH,'a.bed';
while(<FH>){
chomp;
@F=split;
if(exists $normal_chr{$F[0]}){
$fh{$F[0]}->print( "$_\n" );
}else{
print other $_;
}
}
foreach( 1..22 ){
close $fh{$_};
}
close other;
13.JSON格式資料的格式化
題目
學習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 perl
use strict;
use warnings;
use autodie ':all';
use 5.10.0;
use JSON 2;
my $data = from_json( do { local $/; open my $f, '<', $ARGV[0]; scalar <$f> } );
my @fields = qw( technique factor target principal_investigator submission label category type Developmental-Stage organism key );
say join ',', map "\"$_\"", @fields;
for my $item ( @{$data->{items}} ) {
$item->{key} = $item->{label};
no warnings 'uninitialized';
for my $track ( @{$item->{Tracks}} ) {
$item->{label} = $track;
say join ',', map "\"$_\"", @{$item}{@fields};
}
}
參考結果
完成之後的結果應該是：http://biotrainee.com/jbrowse/JBrowse-1.12.1/sample_data/json/modencode/modencodeMetaData.csv
14.多個探針對應一個基因，取平均值、最大值
題目
編寫指令碼對多個探針對應一個基因，取平均值、最大值。
測試資料
R裡面的包自帶很多晶片表達資料，所以測試資料就在下面的R程式碼裡面
R程式碼示例
# 平均值
BiocInstaller::biocLite('CLL')
BiocInstaller::biocLite('hgu95av2.db')
library('hgu95av2.db')
library(CLL)
data(sCLLex)
sCLLex=sCLLex[,1:8] ## 樣本太多，我就取前面8個
group_list=sCLLex$Disease
exprSet=exprs(sCLLex)
exprSet=as.data.frame(exprSet)
exprSet$probe_id=rownames(exprSet)
head(exprSet)
probe2symbol=toTable(hgu95av2SYMBOL)
dat=merge(exprSet,probe2symbol,by='probe_id')
results=t(sapply(split(dat,dat$symbol),function(x) colMeans(x[,1:(ncol(x)-1)])))
15.把counts矩陣轉換成RPKM矩陣
題目
編寫指令碼將hisat2+htseq-counts之後得到reads的counts矩陣轉換成RPKM矩陣。
測試資料
wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_mouse/CCDS.current.txt
grep -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.txt
R程式碼示例
genes_len=read.table("mm10_ccds_length.txt",stringsAsFactors=F)
head(genes_len)
V1 V2
1 -343C11.2 1139
2 00R_AC107638.2 1060
3 00R_Pgap2 1676
4 0610005C13Rik 7363
5 0610006L08Rik 34995
6 0610007P14Rik 9074
colnames(genes_len)<- c("GeneName","Len")
head(exprSet)
GSM860181_priSG-A GSM860182_SG-A GSM860183_SG-B GSM860184_lepSC
00R_AC107638.2 0 1 0 1
0610005C13Rik 20 22 11 27
0610006L08Rik 0 0 0 2
0610007P14Rik 2075 1785 1269 1430
0610009B22Rik 256 376 300 386
0610009E02Rik 22 22 16 28
exprSet<-exprSet[ rownames(exprSet) %in% genes_len$GeneName ,]
total_count<- colSums(exprSet)
neededGeneLength=genes_len[ match(rownames(exprSet), genes_len$GeneName) ,2]
rpkm <- t(do.call( rbind,lapply(1:length(total_count),function(i){
10^9*exprSet[,i]/neededGeneLength/total_count[i]
}) ))
head(rpkm)
rownames(rpkm)=rownames(exprSet)
colnames(rpkm)=colnames(exprSet)
write.table(rpkm,file="rpkm.txt",sep="\t",quote=F)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb=TxDb.Mmusculus.UCSC.mm10.knownGene
gt=transcriptsBy(txdb,by="gene")
lapply(gt[1:40],function(x) max(width(x)))
library(org.Mm.eg.db)
mykeys=
columns(txdb);keytypes(txdb)
neededcols <- c("GENEID", "TXSTRAND", "TXCHROM")
select(txdb, keys=mykeys, columns=neededcols, keytype="TXNAME")
參考結果
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())
## 50 patients and 200 genes
dat=matrix(rnorm(10000,mean=8,sd=4),nrow = 200)
rownames(dat)=paste0('gene_',1:nrow(dat))
colnames(dat)=paste0('sample_',1:ncol(dat))
os_years=abs(floor(rnorm(ncol(dat),mean = 50,sd=20)))
os_status=sample(rep(c('LIVING','DECEASED'),25))
library(survival)
my.surv <- Surv( os_years,os_status=='DECEASED')
## The status indicator, normally 0=alive, 1=dead. Other choices are TRUE/FALSE (TRUE = death) or 1/2 (2=death).
## And most of the time we just care about the time od DECEASED .
fit.KM=survfit(my.surv~1)
fit.KM
plot(fit.KM)
log_rank_p <- apply(dat, 1, function(values1){
group=ifelse(values1>median(values1),'high','low')
kmfit2 <- survfit(my.surv~group)
#plot(kmfit2)
data.survdiff=survdiff(my.surv~group)
p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
})
names(log_rank_p[log_rank_p<0.05])
gender= ifelse(rnorm(ncol(dat))>1,'male','female')
age=abs(floor(rnorm(ncol(dat),mean = 50,sd=20)))
## gender and age are similar with group(by gene expression)
cox_results <- apply(dat , 1, function(values1){
group=ifelse(values1>median(values1),'high','low')
survival_dat <- data.frame(group=group,gender=gender,age=age,stringsAsFactors = F)
m=coxph(my.surv ~ age + gender + group, data = survival_dat)
beta <- coef(m)
se <- sqrt(diag(vcov(m)))
HR <- exp(beta)
HRse <- HR * se
#summary(m)
tmp <- round(cbind(coef = beta, se = se, z = beta/se, p = 1 - pchisq((beta/se)^2, 1),
HR = HR, HRse = HRse,
HRz = (HR - 1) / HRse, HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),
HRCILL = exp(beta - qnorm(.975, 0, 1) * se),
HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)
return(tmp['grouplow',])
})
cox_results[,cox_results[4,]<0.05]
PS: 裡面的os_years應該修改為os_month，因為總的生存期不應該長達幾十年，而且因為ag和os_years都是隨機生成的，可能會出現很不符合自然科學的現象。但是這個是模擬資料，請大家不用較真。
17.對多個差異分析結果直接取交集並集
題目
編寫指令碼對每兩個差異分析結果計算基因交集個數與基因並集個數的比值，得到一個比值矩陣。
測試資料
用perl單行命令模擬資料：
perl -e 'BEGIN{use List::Util qw/shuffle/;@gene=A..Z}{foreach(1..10){@[email protected][(shuffle(0..25))[0..int(rand(20))+4]];print "comparison_$_ <- ";print join(";",@this_genes);print "\n"}}'
總共10次差異分析，每次差異分析得到的差異基因在同一行的後面用大小字母表示。
comparison_1 -> I;G;E;V;C;K;B;W
comparison_2 -> G;E;N;H;Y;M;L;S;K;A;J;O;D;P;R;U;Q;F;Z;C
comparison_3 -> Y;V;U;N;H;K;I;P;S;F;D;X;G;C;Z;J;Q;T;W;O;E;M
comparison_4 -> N;T;K;B;H;Z;W;C;Q;I;V;F;D;S;R;Y;J;X;P;O;G;L;A
comparison_5 -> G;J;A;H;W;T;Z;E;Y;S;R
comparison_6 -> Z;M;D;R;P;G;L;W;Y;U;X;E;A;S;T;I;H
comparison_7 -> H;Z;T;O;W;Q;M;X;J;N;U;K;F;P;I;C;S;Y;A;B
comparison_8 -> A;R;L;T;W;Q;S;F;P;X;E;V;Y;G;K;J;Z;C
comparison_9 -> J;X;K;D;O;H;L;F;C;P;R;N
comparison_10 -> G;S;K;H;C;O;W;F;Q;X
R程式碼示例
a=readLines('test.txt')
n=unlist(lapply(a , function(x){
trimws(strsplit(x,'<-')[[1]][1])
}))
v=lapply(a , function(x){
trimws(strsplit(trimws(strsplit(x,'<-')[[1]][2]),';')[[1]])
})
names(v)=n
tmp=unlist(lapply(v, function(x){
lapply(v, function(y){
p1=length(intersect(x,y))
p2=length(union(x,y))
return(p1/p2)
})
}))
out=matrix(tmp,nrow = length(n))
rownames(out)=n
colnames(out)=n
out[1:5,1:5]
參考結果
結果的前五行如下：
comparison_1 comparison_2 comparison_3 comparison_4 comparison_5
comparison_1 1.0000000 0.1666667 0.3043478 0.2916667 0.1875000
comparison_2 0.1666667 1.0000000 0.6800000 0.6538462 0.4090909
comparison_3 0.3043478 0.6800000 1.0000000 0.7307692 0.3750000
comparison_4 0.2916667 0.6538462 0.7307692 1.0000000 0.4166667
comparison_5 0.1875000 0.4090909 0.3750000 0.4166667 1.0000000
18.根據GTF格式的基因註釋檔案得到人所有基因的染色體座標
題目
從gencode資料庫裡面可以下載所有的gtf檔案，編寫指令碼得到基因的染色體、起始終止座標如下：
[[email protected]]$ head protein_coding.hg19.position
chr1 69091 70008 OR4F5
chr1 367640 368634 OR4F29
chr1 621096 622034 OR4F16
chr1 859308 879961 SAMD11
chr1 879584 894689 NOC2L
chr1 895967 901095 KLHL17
chr1 901877 911245 PLEKHN1
chr1 910584 917473 PERM1
chr1 934342 935552 HES4
chr1 936518 949921 ISG15
[[email protected]]$ head protein_coding.hg38.position
chr1 69091 70008 OR4F5
chr1 182393 184158 FO538757.2
chr1 184923 200322 FO538757.1
chr1 450740 451678 OR4F29
chr1 685716 686654 OR4F16
chr1 923928 944581 SAMD11
chr1 944204 959309 NOC2L
chr1 960587 965715 KLHL17
chr1 966497 975865 PLEKHN1
chr1 975204 982093 PERM1