问题：stringtie组装转录本后会改变原本参考基因组注释文件中的gene_id
-
stringtie 组装、merge的stringtie_merged.gtf注释文件中的gene id 与ballgown生成的gene id 相同
-
stringtie 组装、merge的stringtie_merged.gtf注释文件中的gene id 与prepDE.py生成的counts矩阵的gene id 相同
-
stringtie会将部分参考基因组gtf注释文件的gene_id重命名为MSTRG，并在组装的注释文件中保留为ref_gene_id
思路：
如果想恢复ballgown中的gene_id 或counts矩阵中的gene_id为ref_gene_id，将stringtie merge的注释文件中的gene_id、ref_gene_id提出，在R中使用merge函数进行合并，将ref_gene_id 中的空白单元格赋值为NA（空白单元格即为stringtie新组装的转录本基因），再用id列中的MSTRG基因名将ref_gene_id中的NA替换。
实现：
1）awk提取stringtie组装的stringtie_merged.gtf注释文件中的gene_id和ref_gene_id，sort将首列排序，并保留首列gene_id去重复之后的行。
cat stringtie_merged.gtf | grep "gene_id"| awk '{if($3=="transcript")print}' |awk -F "\t" '{print $9}' | awk -F ";" '{print$1,$3}' |awk '{print $2"\t"$4}'|sed 's/\"//g' | sort -k 1,1 -rk 2,2 -V | sort -uk1,1 -V | less > gene_id_ref_gene_id.txt
还可以加上染色体位置
cat stringtie_merged.gtf | awk '{if($3=="transcript")print}' | awk '{print $1"\t"$4"\t"$5}' | less > gene_position.txt
cat stringtie_merged.gtf | awk '{if($3=="transcript")print}' |awk -F "\t" '{print $9}' | awk -F ";" '{print$1,$3}' |awk '{print $2"\t"$4}'|sed 's/\"//g' | less > gene_id.txt
paste -d "\t" gene_position.txt gene_id.txt | awk -F "\t" '{print $4"\t"$5"\t"$1"\t"$2"\t"$3}' | sort -t $'\t' -k1,1 -rk2,2 -V | sort -t $'\t' -uk1,1 -V | less > gene_id_ref_gene_id.txt
2）R合并results_gene的dataframe和gene_id_ref_gene_id.txt
results_gene = stattest(bg, feature="gene",covariate="condition", getFC=TRUE, meas="FPKM")
genes <- read.table("gene_id_ref_gene_id.txt",header = F,sep = "\t")
names(genes) = c("gene_id","ref_gene_id")
results_gene_1 <- merge(results_gene,genes,by.x = "id",by.y = "gene_id",all = TRUE )
results_gene_1$ref_gene_id[results_gene_1$ref_gene_id == ""] <- "NA"
library("tidyverse")
results_gene_1 <- results_gene_1 %>%
mutate(ref_gene_id = if_else(ref_gene_id == "NA", id, ref_gene_id))
results_gene_1 <- results_gene_1 %>%
select(feature,ref_gene_id,fc,pval,qval)
names(results_gene_1) = c("feature","id","fc","pval","qval")
网友评论