搜索
查看: 3359|回复: 5

TCGA计算FPKM时计算基因长度

[复制链接]

4

主题

11

帖子

85

积分

注册会员

Rank: 2

积分
85
发表于 2018-4-7 19:37:49 | 显示全部楼层 |阅读模式
对于RNA数据TCGA官网提供了三种格式的文件,分别为:
.htseq.counts.gz
.FPKM.txt.gz
.FPKM-UQ.txt.gz
三者之间的关系如下图:

链接地址为:
https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/
从上图我们可以看到从htseq.counts.gz转化到FPKM.txt.gz需要知道映射到所有蛋白质编码基因的读长总和RCpc和基因的长度L。本次我们就主要来说说如何利用官方提供的注释文件计算基因的长度。


TCGA利用的是GRCh38参考基因组,TCGA提供了改注释文件的下载:https://gdc.cancer.gov/about-data/data-harmonization-and-generation/gdc-reference-files
(也可从GENCODE官网下载相应版本)
下载后为gencode.v22.annotation.gtf.gz

该文件具体格式说明见GENCODE官网:
https://www.gencodegenes.org/data_format.html

我们利用awk命令提取了以下几列信息,然后用R语言计算基因长度
"Chromosome_name","Annotiation_source","Feature_type","Genomic_start_location","Genomic_end_location","Gene_ID","Gene_type","Gene_status","Gene_name"

#首先建立一个.awk文件
vim Get_Simple_Annotation.awk
[Shell] 纯文本查看 复制代码
#!/usr/bin/awk -f


[mw_shl_code=perl,true]#!/usr/bin/awk -f

BEGIN{
	print "Chromosome_name","Annotiation_source","Feature_type","Genomic_start_location","Genomic_end_location","Gene_ID","Gene_type","Gene_status","Gene_name" > "Simple_Annotiation.txt"
}

/gene_id/{
	#染色体名字
	chr_name = $1;
	#注释来源
	annot_source = $2;
	#特征类型
	feature_type = $3;
	#基因起始位点
	genomic_start = $4;
	#基因终止位点
	genomic_end = $5;
	#gene_id
	split($9, Nine_Column, ";");
	split(Nine_Column[1], ENGS_ID, "\"");
	Gene_id = ENGS_ID[2];

	#Gene_status
	for (i in Nine_Column)
	{
		if (Nine_Column[i] ~ /^ gene_status/)
		{
			split(Nine_Column[i], gene_status, "\"");
			Gene_status = gene_status[2];
		}
	}

	#Gene_type
	for (i in Nine_Column)
	{
		if(Nine_Column[i] ~ /^ gene_type/)
		{
			split(Nine_Column[i], gene_type, "\"");
			Gene_type = gene_type[2];
		}
	}

	#Gene_name
	for (i in Nine_Column)
	{
		if(Nine_Column[i] ~ /^ gene_name/)
		{
			split(Nine_Column[i], gene_name, "\"");
			Gene_name = gene_name[2];
		}
	}

	#输出结果
	print chr_name,annot_source,feature_type,genomic_start,genomic_end,Gene_id,Gene_type,Gene_status,Gene_name >> "Simple_Annotiation.txt";
}


#输出结果


#利用R语言计算基因长度
#每个基因的长度为所有外显子长度之和,要除去外显子之间的重叠部分,可以把这个问题看作是求多个区间的并集。
[Plain Text] 纯文本查看 复制代码
options(stringsAsFactors = FALSE)

FilePath <- "D:/MyProject/TCGA-CESC/AnnotationFile/Simple_Annotiation.txt"

Original_Annotiation <- read.table(file = FilePath, header = TRUE)

ENGS_ID_Identify <- Original_Annotiation[!duplicated(Original_Annotiation[, "Gene_ID"]), "Gene_ID"]

#创建数据框保存变量结果
Gene_length <- data.frame()

#一个一个处理
for(i in 1 : length(ENGS_ID_Identify)){
  Index <- intersect(which(Original_Annotiation[, "Feature_type"] == "exon"), which(Original_Annotiation[, "Gene_ID"] == ENGS_ID_Identify[i]))
  Target_data <- Original_Annotiation[Index, c("Genomic_start_location", "Genomic_end_location")]
  Target_sort <- Target_data[order(Target_data[, "Genomic_start_location"]), ]
  
  Result_data <- data.frame()
  
  if (nrow(Target_sort) <= 1){
    Result_data[1,1] <- Target_sort[1,1]
    Result_data[1,2] <- Target_sort[1,2]
  }else{
    Result_data[1,1] <- Target_sort[1,1]
    Result_data[1,2] <- Target_sort[1,2]
    for (j in 2 : nrow(Target_sort)) {
      if(Result_data[nrow(Result_data), 2] < Target_sort[j, 1]){
        Result_data[nrow(Result_data) + 1, ] <- Target_sort[j, ]
      }else{
        Result_data[nrow(Result_data), 2] <- max(Target_sort[j, 2], Result_data[nrow(Result_data), 2])
      }
    }
  }
  
  #计算基因长度
  # Length <- 0
  # for (k in 1 : nrow(Result_data)){
  #   Length <- Result_data[k, 2] - Result_data[k, 1] + 1 + Length
  # }
  Length <- sum(Result_data[, 2] - Result_data[, 1]) + nrow(Result_data)
  
  Gene_length[i, 1] <- ENGS_ID_Identify[i]
  Gene_length[i, 2] <- Length
}
colnames(Gene_length) <- c("ENGD_ID", "Gene_length")
write.table(Gene_length, file = "D:/MyProject/TCGA-CESC/AnnotationFile/Gene_length_union.txt", row.names = FALSE, quote = FALSE)

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2018-4-9 20:40:10 | 显示全部楼层
详细!FPKM结果跟官方的一致吗
回复 支持 反对

使用道具 举报

365

主题

512

帖子

1713

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
1713
发表于 2018-4-10 09:27:37 | 显示全部楼层
第一次看到有人写awk脚本呀!!!
厉害
回复 支持 反对

使用道具 举报

4

主题

11

帖子

85

积分

注册会员

Rank: 2

积分
85
 楼主| 发表于 2018-4-10 14:10:43 | 显示全部楼层
anlan 发表于 2018-4-9 20:40
详细!FPKM结果跟官方的一致吗

我比对了一下是一致的,官方从Counts转换到FPKM基因长度不是按照所有外显子长度之和来计算的,而是考虑到外显子之间的重叠。
回复 支持 反对

使用道具 举报

4

主题

47

帖子

768

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
768
发表于 2018-4-22 23:00:40 | 显示全部楼层
这个awk脚本效率有点低。
回复 支持 反对

使用道具 举报

4

主题

47

帖子

768

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
768
发表于 2018-4-23 10:34:28 | 显示全部楼层
对于这种尽量避免使用for循环,需要的时间太长了,我不知道你这个脚本,完成整个过程需要多长时间,但肯定效率不高。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-8-24 22:34 , Processed in 0.031726 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.