問題：stringtie組裝轉(zhuǎn)錄本后會改變原本參考基因組注釋文件中的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
思路：
如果想恢復(fù)ballgown中的gene_id 或counts矩陣中的gene_id為ref_gene_id，將stringtie merge的注釋文件中的gene_id印机、ref_gene_id提出，在R中使用merge函數(shù)進行合并射赛，將ref_gene_id 中的空白單元格賦值為NA（空白單元格即為stringtie新組裝的轉(zhuǎn)錄本基因）多柑，再用id列中的MSTRG基因名將ref_gene_id中的NA替換。
實現(xiàn)：
1）awk提取stringtie組裝的stringtie_merged.gtf注釋文件中的gene_id和ref_gene_id咒劲，sort將首列排序顷蟆，并保留首列g(shù)ene_id去重復(fù)之后的行腐魂。
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/\"http://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/\"http://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")
