听说你需要编程练习题？

首先上目录：

生信编程很简单↩
人类基因组的外显子区域的长度↩
hg19基因组序列的一些探究↩
hg38每条染色体的基因、转录本分布↩
多个同样行列式文件的合并↩
根据GTF画基因的多个转录本结构↩
下载最新版的KEGG信息，并且解析好↩
写超几何分布检验↩
ID转换↩
根据指定染色体及坐标得到序列↩
根据指定染色体及坐标得到位置信息↩
把文件内容按照染色体分开写出↩
JSON格式数据的格式化↩
多个探针对应一个基因，取平均值、最大值↩
把counts矩阵转换成RPKM矩阵↩
对有临床信息的表达矩阵批量做生存分析↩
对多个差异分析结果直接取交集并集↩
根据GTF格式的基因注释文件得到人所有基因的染色体坐标↩

列出我们板块的人气最旺的20个题目，作为第五章节学习后的复习题目。

生信编程很简单[1]

[1]: 生信编程很简单

编程语言系统入门

题目

对FASTQ的操作

• 5,3段截掉几个碱基

• 序列长度分布统计

• FASTQ 转换成 FASTA

• 统计碱基个数及GC%

对FASTA的操作

• 取互补序列

• 取反向序列

• DNA to RNA

• 大小写字母形式输出

• 每行指定长度输出序列

• 按照序列长度/名字排序

• 提取指定ID的序列

• 随机抽取序列

高级难度

• 根据坐标取序列

• 多文件合并

• 根据ID列表取序列

• GTF文件探索

• 简并碱基的引物序列还原成多条序列

• snp进行注释并格式化输出

下载安装bowtie2（内含测试数据）

1. cd ~/biosoft

2. mkdir bowtie &&  cd bowtie

3. wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.2.9/bowtie2-2.2.9-linux-x86_64.zip  

4. unzip bowtie2-2.2.9-linux-x86_64.zip

人类基因组的外显子区域的长度[^2]

[2]: 人类基因组的外显子区域的长度

题目

下载人类外显子的坐标文件，编写代码统计外显子区域的长度。

测试数据

• Rbioconductor的TxDb.Hsapiens.UCSC.hg19.knownGene包

• NCBI数据库：ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_human/

R实现代码示例

1. rm(list=ls())

2. a=read.table(choose.files(),sep = '\t',stringsAsFactors = F,header = T) # 选择你下的CCDs文件

3. tmp <- apply(a[1:100,], 1, function(gene){ # 取前100行数据分析调试

4. # gene=a[3,]

5.  x=gene[10] # Column10 外显子坐标位置列

6.  if(grepl('\\]',x)){ # 判断x中是否存在有]这样的符号，如果有就利用正则替换掉。

7.    x=sub('\\[','',x)

8.    x=sub('\\]','',x)

9.    # 这个时候得到的对象还是像这样的“880073-880179, 880436-880525……”

10.    tmp <- strsplit(as.character(x),',')[[1]]# 我们先从逗号开始分割成小块

11.    start <- as.numeric(unlist(lapply(tmp,function(y){# 取开始位点

12.      strsplit(as.character(y),'-')[[1]][1]

13.    })))

14.    end <- as.numeric(unlist(lapply(tmp,function(y){ # 取结束位点

15.      strsplit(as.character(y),'-')[[1]][2]

16.    })))

17.    gene_d <- data.frame(gene=gene[3], # 将基因名，染色体，开始、结束位点绑定为数据框

18.                      chr=gene[1],

19.                      start=start,

20.                      end=end

21.    )

22.    return (gene_d)#返回数据框

23.  }

24. })

25. tmp_pos=c() # 构造一个空的向量

26. lapply(tmp[1:10], function(x){ # 取前10个list文件计算调试

27. # print(x)

28.  if ( !is.null(x)){

29.    apply(x, 1,function(y){

30.      #print(y)

31.      for ( i in as.numeric(y[3]):as.numeric(y[4]) ) # y[3]为坐标起点，y[4]为终止坐标，历编

32.        tmp_pos <<- c(tmp_pos,paste0(y[2],"-",i))

33.    })

34.  }

35. })

36. length(tmp_pos) # 计算exon的长度

37. length(unique(tmp_pos)) # 计算去重后的exon的长度

hg19基因组序列的一些探究[^3]

[3]: hg19基因组序列的一些探究

题目

求：hg19 每条染色体长度，每条染色体N的含量，GC含量。

（高级作业：蛋白编码区域的GC含量会比基因组其它区域的高吗？ ）

测试数据

• hg19基因组序列下载

1. wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz # 也可以在浏览器上下载

2. tar xvzf chromFa.tar.gz

3. cat *.fa > hg19.fa

4. rm chr*.fa # 先把着急删，我待会可能要那他测试运行速度

• 简单测试数据

1. >chr_1

2. ATCGTCGaaAATGAANccNNttGTA

3. AGGTCTNAAccAAttGggG

4. >chr_2

5. ATCGAATGATCGANNNGccTA

6. AGGTCTNAAAAGG

7. >chr_3

8. ATCGTCGANNNGTAATggGA

9. AGGTCTNAAAAGG

10. >chr_4

11. ATCGTCaaaGANNAATGANGgggTA

Perl代码示例

单行命令 {-}

1. 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

完整代码 {-}

1. while (<>){

2. chomp;

3. if(/^>/){

4. $chr=$_

5. }

6. else{

7. $A_count{$chr}+=($_=~tr/Aa//);

8. $T_count{$chr}+=($_=~tr/Tt//);

9. $C_count{$chr}+=($_=~tr/Cc//);

10. $G_count{$chr}+=($_=~tr/Gg//);

11. $N_count{$chr}+=($_=~tr/Nn//);

12. }

13. }

14. foreach (sort keys  %N_count){

15. $length = $A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_}+$N_count{$_};

16. $gc = ($G_count{$_}+$C_count{$_})/($A_count{$_}+$T_count{$_}+$C_count{$_}+$G_count{$_});

17. print "$_\t$A_count{$_}\t$T_count{$_}\t$C_count{$_}\t$G_count{$_}\t$N_count{$_}\t$length\t$gc\n"

18. }

参考结果{-}

结果如下；

1. >chr_1        13        10        7        10        4

2. >chr_2        11        6        5        8        4

3. >chr_3        10        6        3        10        4

4. >chr_4        9        4        2        7        3

hg38每条染色体的基因、转录本分布[^4]

[4]: hg38每条染色体的基因、转录本分布

题目

对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的基因？它们的特征又有哪些变化呢？

测试数据

1. wget -c ftp://ftp.ensembl.org/pub/release-87/gtf/homo_sapiens/Homo_sapiens.GRCh38.87.chr.gtf.gz

2. gzip -d Homo_sapiens.GRCh38.87.chr.gtf.gz

代码示例

1. # 每条染色体的基因个数

2. zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{print if $F[2] eq "gene" }' |cut -f 1 |sort |uniq -c

3. # 基因分类

4. zcat Homo_sapiens.GRCh38.87.chr.gtf.gz |perl -alne '{next unless $F[2] eq "gene" ;/gene_biotype "(.*?)";/;print $1}'  |sort |uniq -c

多个同样行列式文件的合并[^5]

[5]: 多个同样行列式文件的合并

题目

将htseq-count生成的所有独立样本文件进行合并（每个样品对应一个文件，包括了所有基因表达量）。 希望通过编程处理每个文件得到输出的表达矩阵（行名是基因名，列名是样品名），如下所示：

1

测试数据

模拟数据

用perl脚本模仿htseq-count计算每个样本所有的基因表达量的输出独立样本文件：

1. ## 首先新建文件tmp.sh 输入这个代码：

2. perl -le '{print "gene_$_\t".int(rand(1000)) foreach 1..99}'

3. ## 然后用perl脚本调用这个tmp.sh文件：

4. perl -e 'system(" bash tmp.sh >$_.txt") foreach a..z'

5. ## 这样就生成了a~z这26个样本的counts文件

第一列是基因，第二列是该基因的counts值，共有a~z这26个样本的counts文件，需要合并成一个大的行列式，这样才能导入到R里面做差异分析，如果手工用excel表格做，当然是可以的，但是太麻烦，如果有500个样本，正常人都不会去手工做了，需要编程。

每个样本的基因顺序并不一致，这时候你应该怎么做呢？

真实数据

实际需求如下：GSE48213里面有56个文件，需要合并成一个表达矩阵，来根据cell-line的不同，分组做差异分析。可以查看paper

1. wget -c ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE48nnn/GSE48213/suppl/GSE48213_RAW.tar

2. tar -xf GSE48213_RAW.tar

3. gzip -d *.gz

代码示例

1. ## 首先在GSE48213_RAW目录里面生成tmp.txt文件，使用shell脚本：

2. awk '{print FILENAME"\t"$0}' * |grep -v EnsEMBL_Gene_ID >tmp.txt

3. ## 然后把tmp.txt导入R语言里面用reshape2处理即可！

4. setwd('tmp/GSE48213_RAW/')

5. a=read.table('tmp.txt',sep = '\t',stringsAsFactors = F)

6. library(reshape2)

7. fpkm <- dcast(a,formula = V2~V1)

根据GTF画基因的多个转录本结构[^6]

[6]: 根据GTF画基因的多个转录本结构

题目

从NCBI,ENSEMBL,UCSC,GENCODE数据库下载各种GTF注释文件，编写代码得到所有基因的转录本个数，以及每个转录本的外显子的坐标，绘制如下转录本结构图：

2

比如对这个ANXA1基因来说，非常多的转录本，但是基因的起始终止坐标，是所有转录本起始终止坐标的极大值和极小值。同时，它是一个闭合基因，因为它存在一个转录本的起始终止坐标等于该基因的起始终止坐标。可以看到它的外显子并不是非常整齐的，虽然多个转录本会共享某些外显子，但是也存在某些外显子比同区域其它外显子长的现象。

测试数据

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

2. gzip -d gencode.v7.annotation_goodContig.gtf.gz

R实现代码示例

1. rm(list=ls())

2. ## http://www.broadinstitute.org/cancer/cga/sites/default/files/data/tools/rnaseqc/gencode.v7.annotation_goodContig.gtf.gz

3. setwd('tmp')

4. gtf <- read.table('gencode.v7.annotation_goodContig.gtf.gz',stringsAsFactors = F,

5.                  header = F,comment.char = "#",sep = '\t'

6.                  )

7. table(gtf[,2])

8. gtf <- gtf[gtf[,2] =='HAVANA',]

9. gtf <- gtf[grepl('protein_coding',gtf[,9]),]

10. lapply(gtf[1:10,9], function(x){

11.  y=strsplit(x,';')

12. })

13. gtf$gene <- lapply(gtf[,9], function(x){

14.  y <- strsplit(x,';')[[1]][5]

15.  strsplit(y,'\\s')[[1]][3]

16.  }

17. )

18. draw_gene = 'TP53'

19. structure = gtf[gtf$gene==draw_gene,]

20. colnames(structure) =c(

21.  'chr','db','record','start','end','tmp1','tmp2','tmp3','tmp4','gene'

22. )

23. gene_start <- min(c(structure$start,structure$end))

24. gene_end <- max(c(structure$start,structure$end))

25. tmp_min=min(c(structure$start,structure$end))

26. structure$new_start=structure$start-tmp_min

27. structure$new_end=structure$end-tmp_min

28. tmp_max=max(c(structure$new_start,structure$new_end))

29. num_transcripts=nrow(structure[structure$record=='transcript',])

30. tmp_color=rainbow(num_transcripts)

31. x=1:tmp_max;y=rep(num_transcripts,length(x))

32. #x=10000:17000;y=rep(num_transcripts,length(x))

33. plot(x,y,type = 'n',xlab='',ylab = '',ylim = c(0,num_transcripts+1))

34. title(main = draw_gene,sub = paste("chr",structure$chr,":",gene_start,"-",gene_end,sep=""))

35. j=0;

36. tmp_legend=c()

37. for (i in 1:nrow(structure)){

38.  tmp=structure[i,]

39.  if(tmp$record == 'transcript'){

40.    j=j+1

41.    tmp_legend=c(tmp_legend,paste("chr",tmp$chr,":",tmp$start,"-",tmp$end,sep=""))

42.  }

43.  if(tmp$record == 'exon') lines(c(tmp$new_start,tmp$new_end),c(j,j),col=tmp_color[j],lwd=4)

44. }

参考结果

3

下载最新版的KEGG信息，并且解析好[^7]

[7]: 下载最新版的KEGG信息，并且解析好

题目

下载最新版的KEGG注释文本文件，编写脚本整理成kegg的pathway的ID与基因ID的对应格式。

测试数据

1 首先打开KEGG官方网站，网页中展示出了各个物种的分类、拉丁名称、英文名称等信息。

4

2 直接网页中搜索（Ctrl + F）需要下载的物种英文名称或拉丁名。如果不确定物种名称，网站中提供了详细的分类系统，也可根据前面的物种分类信息进行查找。 本文以拟南芥为例，搜索“Arabidopsis thaliana”即可找到。找到后点击物种名称前的3个字母缩写链接（下图红色框中的位置）。

3 进入后的网页中包含了物种的一些基因组信息，点击上方的“Brite hierarchy”，进入后再点击“KEGG Orthology (KO)”；

6

4 在跳转出的网页中点击“Download htext”，弹出下载窗口进行下载，国外网站有时会出现无法下载的情况，多试几次即可；

5 当然，下载好之后还没有结束。下载得到文本文件，可以看到里面的结构层次非常清楚，C开头的就是kegg的pathway的ID所在行，D开头的就是属于它的kegg的所有的基因。A,B是kegg的分类，总共是6个大类，42个小类。

• 8

Perl代码示例

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

参考结果

9写超几何分布检验[^8]

[8]: 写超几何分布检验

题目

学习GO/KEGG的富集分析的原理，编写代码实现超几何分布检验，将得到的结果与测试数据中的kegg.enrichment.html进行比较。

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

测试数据

• kegg2gene（第六讲kegg数据解析结果）

暂时不用最新版的kegg注释数据，为了能够统一答案

• 差异基因list和背景基因list(R代码)

1. source("http://bioconductor.org/biocLite.R")

2. biocLite("org.Hs.eg.db")

3. biocLite("KEGG.db")

4. biocLite("GOstats")

5. biocLite("hgu95av2.db")

6. library(org.Hs.eg.db)

7. library(KEGG.db)

8. library(GOstats)

9. library("hgu95av2.db")

10. ##得到kegg2gene.list(KEGG注释信息)

11. tmp=toTable(org.Hs.egPATH)

12. write.table(tmp,'kegg2gene.list.txt',quote = F,row.names = F)

13. ## 得到universeGeneIds.txt(背景基因list)

14. ls('package:hgu95av2.db')

15. universeGeneIds <- unique(mappedRkeys(hgu95av2ENTREZID))

16. write.table(universeGeneIds,'universeGeneIds.txt',quote = F,row.names = F)

17. ## 得到diff_gene.txt(差异基因list)

18. set.seed('123456.789')

19. diff_gene = sample(universeGeneIds,300)

20. write.table(diff_gene,'diff_gene.txt',quote = F,row.names = F)

21. ## 得到kegg.enrichment.html(富集分析结果)

22. annotationPKG='org.Hs.eg.db'

23. hyperG.params = new("KEGGHyperGParams", geneIds=diff_gene, universeGeneIds=universeGeneIds, annotation=annotationPKG,

24. categoryName="KEGG", pvalueCutoff=1, testDirection = "over")

25. KEGG.hyperG.results = hyperGTest(hyperG.params);

26. htmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))

10

R代码示例

下载链接： https://github.com/jmzeng1314/humanid/blob/master/R/hyperGtest_jimmy.R

1. library("hgu95av2.db")[/align][align=left]ls('package:hgu95av2.db')

2. universeGeneIds <- unique(mappedRkeys(hgu95av2ENTREZID))

3. set.seed('123456.789')

4. diff_gene = sample(universeGeneIds,300)

5. library(org.Hs.eg.db)

6. library(KEGG.db)

7. tmp=toTable(org.Hs.egPATH)

8. GeneID2kegg<<- tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)

9. kegg2GeneID<<- tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)

10. #results <- hyperGtest_jimmy(GeneID2kegg,kegg2GeneID,diff_gene,universeGeneIds)

11. GeneID2Path=GeneID2kegg

12. Path2GeneID=kegg2GeneID

13. diff_gene_has_path=intersect(diff_gene,names(GeneID2Path))

14. n=length(diff_gene) #306

15. N=length(universeGeneIds) #5870

16. results=c()

17. for (i in names(Path2GeneID)){

18.  #i="04672"

19.  M=length(intersect(Path2GeneID[[i]],universeGeneIds))

20.  #print(M)

21.  if(M<5)

22.    next

23.  exp_count=n*M/N

24.  #print(paste(n,N,M,sep="\t"))

25.  k=0

26.  for (j in diff_gene_has_path){

27.    if (i %in% GeneID2Path[[j]]) k=k+1

28.  }

29.  OddsRatio=k/exp_count

30.  if (k==0) next

31.  p=phyper(k-1,M, N-M, n, lower.tail=F)

32.  #print(paste(i,p,OddsRatio,exp_count,k,M,sep="    "))

33.  results=rbind(results,c(i,p,OddsRatio,exp_count,k,M))

34. }

35. colnames(results)=c("PathwayID","Pvalue","OddsRatio","ExpCount","Count","Size")

36. results=as.data.frame(results,stringsAsFactors = F)

37. results$p.adjust = p.adjust(results$Pvalue,method = 'BH')

38. results=results[order(results$Pvalue),]

39. rownames(results)=1:nrow(results)

ID转换[^9]

[9]: ID转换

题目

probeid，geneid，gene_name, symbol之间的转换。

测试数据

• 需要转换的探针ID(变量my_probe)

1. rm(list=ls())

2. library("hgu95av2.db")

3. ls('package:hgu95av2.db')

4. probe2entrezID=toTable(hgu95av2ENTREZID)

5. probe2symbol=toTable(hgu95av2SYMBOL)

6. probe2genename=toTable(hgu95av2GENENAME)

7. my_probe = sample(unique(mappedLkeys(hgu95av2ENTREZID)),30)

8. tmp1 = probe2symbol[match(my_probe,probe2symbol$probe_id),]

9. tmp2 = probe2entrezID[match(my_probe,probe2entrezID$probe_id),]

10. tmp3 = probe2genename[match(my_probe,probe2genename$probe_id),]

11. write.table(my_probe,'my_probe.txt',quote = F,col.names = F,row.names =F)

12. write.table(tmp1$symbol,'my_symbol.txt',quote = F,col.names = F,row.names =F)

13. 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文件。

1. library("illuminaHumanv4.db")

2. ls('package:illuminaHumanv4.db')

3. probe2entrezID=toTable(illuminaHumanv4ENTREZID)

4. probe2symbol=toTable(illuminaHumanv4SYMBOL)

5. probe2genename=toTable(illuminaHumanv4GENENAME)

6. my_probe = sample(unique(mappedLkeys(illuminaHumanv4ENTREZID)),30)

7. probe2symbol[match(my_probe,probe2symbol$probe_id),]

8. probe2entrezID[match(my_probe,probe2entrezID$probe_id),]

9. probe2genename[match(my_probe,probe2genename$probe_id),]

R代码示例

基因的转换：运行下面的R代码，得到的mysymbolgene和myentrezgene就是需要转换的ID。

1. library("illuminaHumanv4.db")

2. ls('package:illuminaHumanv4.db')

3. my_entrez_gene = sample(unique(mappedRkeys(illuminaHumanv4ENTREZID)),30)

4. my_symbol_gene = sample(unique(mappedRkeys(illuminaHumanv4SYMBOL)),30)

5. library("org.Hs.eg.db")

6. ls('package:org.Hs.eg.db')

7. entrezID2symbol <- toTable(org.Hs.egSYMBOL)

8. entrezID2symbol[match(my_entrez_gene,entrezID2symbol$gene_id),]

9. entrezID2symbol[match(my_symbol_gene,entrezID2symbol$symbol),]

根据指定染色体及坐标得到序列[^10]

[10]: 根据指定染色体及坐标得到序列

题目

参考基因组hg19，指定染色体及坐标"chr5","8397384"，编写程序得到这个坐标以及它上下一个碱基。（机器无法计算hg19，则使用测试数据，指定坐标是 3号染色体的第6个碱基。）

测试数据

1. >chr_1

2. ATCGTCGaaAATGAANccNNttGTA

3. AGGTCTNAAccAAttGggG

4. >chr_2

5. ATCGAATGATCGANNNGccTA

6. AGGTCTNAAAAGG

7. >chr_3

8. ATCGTCGANNNGTAATggGA

9. AGGTCTNAAAAGG

10. >chr_4

11. ATCGTCaaaGANNAATGANGgggTA

根据指定染色体及坐标得到位置信息[^11]

[11]: 根据指定染色体及坐标得到位置信息

题目

基因的chr,start,end都是已知的（坐标是hg38系统），任意给定基因组的chr:pos(chr1:2075000-2930999), 判断该区间在哪个基因上面？（可用现成软件bedtools）

测试数据

1. chr7        148697841        148698941

2. chr7        148698942        148699029

3. chr7        148699911        148701053

4. chr7        148701109        148701307

5. chr7        148701354        148702694

6. chr7        148703100        148703520

7. chr7        148703831        148704175

8. chr7        148704484        148704734

9. chr7        148704857        148705937

10. chr7        148706271        148706671

把文件内容按照染色体分开写出[^12]

[12]: 把文件内容按照染色体分开写出

题目

根据染色体把文件拆分成1~22和其它染色体的两个文件呢？

测试数据

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

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

1. chr2    43995310    43995986

2. chr17  49788603    49789067

3. chr17  59565573    59566163

4. chr19  8390308 8390745

5. chr12  49188033    49189033

6. chr7    974903  975570

7. chr7    98878532    98879500

8. chr7    44044672    44045322

9. chr1    153634052  153634772

10. chr11  60905850    60906575

Perl代码示例

1. use FileHandle;

2. foreach( 1..22 ){

3.    $normal_chr{"chr".$_}=1 ;

4.    open $fh{"chr".$_},">chr$_.txt" or die;

5. }

6. open other,">chr_other.txt" or die;

7. open FH,'a.bed';

8. while(<FH>){

9.    chomp;

10.    @F=split;

11.    if(exists $normal_chr{$F[0]}){

12.        $fh{$F[0]}->print( "$_\n" );

13.    }else{

14.        print other $_;

15.    }

16. }

17. foreach( 1..22 ){

18.    close $fh{$_};

19. }

20. close other;

JSON格式数据的格式化[^13]

[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代码示例

1. #!/usr/bin/env perl

2. use strict;

3. use warnings;

4. use autodie ':all';

5. use 5.10.0;

6. use JSON 2;

7. my $data = from_json( do { local $/; open my $f, '<', $ARGV[0]; scalar <$f> } );

8. my @fields = qw( technique factor target principal_investigator submission label category type Developmental-Stage organism key );

9. say join ',', map "\"$_\"", @fields;

10. for my $item ( @{$data->{items}} ) {

11.    $item->{key} = $item->{label};

12.    no warnings 'uninitialized';

13.    for my $track ( @{$item->{Tracks}} ) {

14.        $item->{label} = $track;

15.        say join ',', map "\"$_\"", @{$item}{@fields};

16.    }

17. }

参考结果

完成之后的结果应该是：http://biotrainee.com/jbrowse/JBrowse-1.12.1/sample_data/json/modencode/modencodeMetaData.csv11

多个探针对应一个基因，取平均值、最大值 [^14]

[14]: 多个探针对应一个基因，取平均值、最大值

题目

编写脚本对多个探针对应一个基因，取平均值、最大值。

测试数据

R代码示例

1. # 平均值

2. BiocInstaller::biocLite('CLL')

3. BiocInstaller::biocLite('hgu95av2.db')

4. library('hgu95av2.db')

5. library(CLL)

6. data(sCLLex)

7. sCLLex=sCLLex[,1:8] ## 样本太多，我就取前面8个

8. group_list=sCLLex$Disease

9. exprSet=exprs(sCLLex)

10. exprSet=as.data.frame(exprSet)

11. exprSet$probe_id=rownames(exprSet)

12. head(exprSet)

13. probe2symbol=toTable(hgu95av2SYMBOL)

14. dat=merge(exprSet,probe2symbol,by='probe_id')

15. results=t(sapply(split(dat,dat$symbol),function(x) colMeans(x[,1:(ncol(x)-1)])))

把counts矩阵转换成RPKM矩阵[^15]

[15]: 把counts矩阵转换成RPKM矩阵

题目

编写脚本将hisat2+htseq-counts之后得到reads的counts矩阵转换成RPKM矩阵。

测试数据

1. wget ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/current_mouse/CCDS.current.txt

2. 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代码示例

1. genes_len=read.table("mm10_ccds_length.txt",stringsAsFactors=F)

2. head(genes_len)

3.              V1    V2

4. 1      -343C11.2  1139

5. 2 00R_AC107638.2  1060

6. 3      00R_Pgap2  1676

7. 4  0610005C13Rik  7363

8. 5  0610006L08Rik 34995

9. 6  0610007P14Rik  9074

10. colnames(genes_len)<- c("GeneName","Len")

11. head(exprSet)

12.              GSM860181_priSG-A GSM860182_SG-A GSM860183_SG-B GSM860184_lepSC

13. 00R_AC107638.2                0              1              0              1

14. 0610005C13Rik                20            22            11              27

15. 0610006L08Rik                  0              0              0              2

16. 0610007P14Rik              2075          1785          1269            1430

17. 0610009B22Rik                256            376            300            386

18. 0610009E02Rik                22            22            16              28

19. exprSet<-exprSet[ rownames(exprSet) %in% genes_len$GeneName ,]

20. total_count<- colSums(exprSet)

21. neededGeneLength=genes_len[  match(rownames(exprSet), genes_len$GeneName) ,2]

22. rpkm <- t(do.call( rbind,lapply(1:length(total_count),function(i){

23.        10^9*exprSet[,i]/neededGeneLength/total_count[i]

24. }) ))

25. head(rpkm)

26. rownames(rpkm)=rownames(exprSet)

27. colnames(rpkm)=colnames(exprSet)

28. write.table(rpkm,file="rpkm.txt",sep="\t",quote=F)

29. library(TxDb.Mmusculus.UCSC.mm10.knownGene)

30. txdb=TxDb.Mmusculus.UCSC.mm10.knownGene

31. gt=transcriptsBy(txdb,by="gene")

32. lapply(gt[1:40],function(x) max(width(x)))

33. library(org.Mm.eg.db)

34. mykeys=

35. columns(txdb);keytypes(txdb)

36. neededcols <- c("GENEID", "TXSTRAND", "TXCHROM")

37. select(txdb, keys=mykeys, columns=neededcols, keytype="TXNAME")

参考结果

对有临床信息的表达矩阵批量做生存分析[^16]

[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等等)的影响。

代码示例

1. rm(list=ls())

2. ## 50 patients and 200 genes

3. dat=matrix(rnorm(10000,mean=8,sd=4),nrow = 200)

4. rownames(dat)=paste0('gene_',1:nrow(dat))

5. colnames(dat)=paste0('sample_',1:ncol(dat))

6. os_years=abs(floor(rnorm(ncol(dat),mean = 50,sd=20)))

7. os_status=sample(rep(c('LIVING','DECEASED'),25))

8. library(survival)

9. my.surv <- Surv( os_years,os_status=='DECEASED')

10. ## The status indicator, normally 0=alive, 1=dead. Other choices are TRUE/FALSE (TRUE = death) or 1/2 (2=death).

11. ## And most of the time we just care about the time od DECEASED .

12. fit.KM=survfit(my.surv~1)

13. fit.KM

14. plot(fit.KM)

15. log_rank_p <- apply(dat, 1, function(values1){

16.  group=ifelse(values1>median(values1),'high','low')

17.  kmfit2 <- survfit(my.surv~group)

18.  #plot(kmfit2)

19.  data.survdiff=survdiff(my.surv~group)

20.  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)

21. })

22. names(log_rank_p[log_rank_p<0.05])

23. gender= ifelse(rnorm(ncol(dat))>1,'male','female')

24. age=abs(floor(rnorm(ncol(dat),mean = 50,sd=20)))

25. ## gender and age are similar with group(by gene expression)

26. cox_results <- apply(dat , 1, function(values1){

27.  group=ifelse(values1>median(values1),'high','low')

28.  survival_dat <- data.frame(group=group,gender=gender,age=age,stringsAsFactors = F)

29.  m=coxph(my.surv ~ age + gender + group, data =  survival_dat)

30.  beta <- coef(m)

31.  se <- sqrt(diag(vcov(m)))

32.  HR <- exp(beta)

33.  HRse <- HR * se

34.  #summary(m)

35.  tmp <- round(cbind(coef = beta, se = se, z = beta/se, p = 1 - pchisq((beta/se)^2, 1),

36.                    HR = HR, HRse = HRse,

37.                    HRz = (HR - 1) / HRse, HRp = 1 - pchisq(((HR - 1)/HRse)^2, 1),

38.                    HRCILL = exp(beta - qnorm(.975, 0, 1) * se),

39.                    HRCIUL = exp(beta + qnorm(.975, 0, 1) * se)), 3)

40.  return(tmp['grouplow',])

41. })

42. cox_results[,cox_results[4,]<0.05]

PS: 里面的osyears应该修改为osmonth，因为总的生存期不应该长达几十年，而且因为ag和os_years都是随机生成的，可能会出现很不符合自然科学的现象。但是这个是模拟数据，请大家不用较真。

对多个差异分析结果直接取交集并集[^17]

[17]: 对多个差异分析结果直接取交集并集

题目

编写脚本对每两个差异分析结果计算基因交集个数与基因并集个数的比值，得到一个比值矩阵。

测试数据

用perl单行命令模拟数据：

1. 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"}}'

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

1. comparison_1 -> I;G;E;V;C;K;B;W

2. comparison_2 -> G;E;N;H;Y;M;L;S;K;A;J;O;D;P;R;U;Q;F;Z;C

3. comparison_3 -> Y;V;U;N;H;K;I;P;S;F;D;X;G;C;Z;J;Q;T;W;O;E;M

4. 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

5. comparison_5 -> G;J;A;H;W;T;Z;E;Y;S;R

6. comparison_6 -> Z;M;D;R;P;G;L;W;Y;U;X;E;A;S;T;I;H

7. comparison_7 -> H;Z;T;O;W;Q;M;X;J;N;U;K;F;P;I;C;S;Y;A;B

8. comparison_8 -> A;R;L;T;W;Q;S;F;P;X;E;V;Y;G;K;J;Z;C

9. comparison_9 -> J;X;K;D;O;H;L;F;C;P;R;N

10. comparison_10 -> G;S;K;H;C;O;W;F;Q;X

R代码示例

1. a=readLines('test.txt')

2. n=unlist(lapply(a , function(x){

3.  trimws(strsplit(x,'<-')[[1]][1])

4. }))

5. v=lapply(a , function(x){

6.  trimws(strsplit(trimws(strsplit(x,'<-')[[1]][2]),';')[[1]])

7. })

8. names(v)=n

9. tmp=unlist(lapply(v, function(x){

10.  lapply(v, function(y){

11.    p1=length(intersect(x,y))

12.    p2=length(union(x,y))

13.    return(p1/p2)

14.  })

15. }))

16. out=matrix(tmp,nrow = length(n))

17. rownames(out)=n

18. colnames(out)=n

19. out[1:5,1:5]

参考结果

结果的前五行如下：

1. comparison_1 comparison_2 comparison_3 comparison_4 comparison_5

2. comparison_1    1.0000000    0.1666667    0.3043478    0.2916667    0.1875000

3. comparison_2    0.1666667    1.0000000    0.6800000    0.6538462    0.4090909

4. comparison_3    0.3043478    0.6800000    1.0000000    0.7307692    0.3750000

5. comparison_4    0.2916667    0.6538462    0.7307692    1.0000000    0.4166667

6. comparison_5    0.1875000    0.4090909    0.3750000    0.4166667    1.0000000

根据GTF格式的基因注释文件得到人所有基因的染色体坐标[^18]

[18]: 根据GTF格式的基因注释文件得到人所有基因的染色体坐标

题目

从gencode数据库里面可以下载所有的gtf文件，编写脚本得到基因的染色体、起始终止坐标如下：

1. [jianmingzeng@gencode]$ head protein_coding.hg19.position

2. chr1        69091        70008        OR4F5

3. chr1        367640        368634        OR4F29

4. chr1        621096        622034        OR4F16

5. chr1        859308        879961        SAMD11

6. chr1        879584        894689        NOC2L

7. chr1        895967        901095        KLHL17

8. chr1        901877        911245        PLEKHN1

9. chr1        910584        917473        PERM1

10. chr1        934342        935552        HES4

11. chr1        936518        949921        ISG15

12. [jianmingzeng@gencode]$ head protein_coding.hg38.position

13. chr1        69091        70008        OR4F5

14. chr1        182393        184158        FO538757.2

15. chr1        184923        200322        FO538757.1

16. chr1        450740        451678        OR4F29

17. chr1        685716        686654        OR4F16

18. chr1        923928        944581        SAMD11

19. chr1        944204        959309        NOC2L

20. chr1        960587        965715        KLHL17

21. chr1        966497        975865        PLEKHN1

22. chr1        975204        982093        PERM1

猜你喜欢

菜鸟入门

数据分析

编程实践

直播基因组分析

生信技能树