最近做基因的miRNA预测，需要用到基因的3'utr序列，可以使用biomaRt包批量获取基因的3'utr序列。

1.安装R包

BiocManager::install("biomaRt")
BiocManager::install("Biostrings")

2. 加载R包

library(biomaRt)　# 用于获取基因的3'utr序列
library(Biostrings) # 将获取的3'utr序列保存成fasta格式
library(clusterProfiler)
library(org.Hs.eg.db)

3. 将基因集转换成ensembl ID

gene = c("TP53","EGFR","VEGFA","APOE","IL6","TGFB1","MTHFR")
gene_table = bitr(gene,fromType = "SYMBOL",toType = "ENSEMBL",OrgDb = "org.Hs.eg.db")
gene_id = gene_table$ENSEMBL
gene_id

4.使用biomaRt包获取基因的3'utr序列

# 指定Ensembl数据库和基因ID
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# 获取基因的 transcript ID
transcript_info <- getBM(attributes = c("ensembl_transcript_id"),
                         filters = "ensembl_gene_id",
                         values = gene_id,
                         mart = ensembl)

# 将 transcript ID提取出来
transcript_id <- transcript_info$ensembl_transcript_id
transcript_id

# 获取 3'UTR 序列
utr_seq <- getSequence(id = transcript_id,
                       type = "ensembl_transcript_id",
                       seqType = "3utr", # 可以更改其参数获取目的序列
                       mart = ensembl)
show(utr_seq)

seqType参数如下

'cdna': for nucleotide sequences

'peptide': for protein sequences

'3utr': for 3' UTR sequences

'5utr': for 5' UTR sequences

'gene_exon': for exon sequences only

'transcript_exon_intron': gives the full unspliced transcript, that is exons + introns

'gene_exon_intron' gives the exons + introns of a gene;'coding' gives the coding sequence only

'coding_transcript_flank': gives the flanking region of the transcript including the UTRs, this must be accompanied with a given value for the upstream or downstream attribute

'coding_gene_flank': gives the flanking region of the gene including the UTRs, this must be accompanied with a given value for the upstream or downstream attribute

'transcript_flank': gives the flanking region of the transcript exculding the UTRs, this must be accompanied with a given value for the upstream or downstream attribute

'gene_flank': gives the flanking region of the gene excluding the UTRs, this must be accompanied with a given value for the upstream or downstream attribute

保存结果
保存成csv文件，第一列是3'utr序列，第二列是ensembl ID。一个基因有多个转录本，并不是所有的转录本都有3'utr序列，可以通过R语言或者直接通过excel删除这部分没有序列的转录本。

write.csv(as.data.frame(utr_seq), "utr_seq.csv", row.names = F)

5.对结果进一步处理

由于对miRNA预测需要使用基因名和3'utr序列，而不是使用转录本ID，因此，需要进一步处理。
使用biomaRt包将转录本ID转换成基因名，不能使用clusterprofiler包的bitr函数进行转换。

# 读取文件
data = read.csv("utr_seq.csv", header = T)
dim(data) 

# 去掉没有3'utr序列的基因
data = subset(data, X3utr != "Sequence unavailable")
dim(data)

# 把转录本ID转换成基因名
gene_name<-getBM(attributes=c("ensembl_transcript_id","external_gene_name"),
                 filters = "ensembl_transcript_id",
                 values = data$ensembl_transcript_id,
                 mart = ensembl)
head(gene_name)

# 合并数据框
df = merge(data,gene_name, by = "ensembl_transcript_id")
head(df)


# 将数据框转换为FASTA格式
fasta_content <- paste(sprintf(">%s\n%s", df$external_gene_name, df$X3utr), collapse = "\n")
# 保存为FASTA文件
writeLines(fasta_content, "gene_3utr.fasta")