lshepard

我想要得到我拥有的EntrezIDs的基因符号，以及这些基因的启动子区域。然而，我不知道我做的是否正确。以下是我所做的:

library(AnnotationHub) library("org.Hs.eg.db") hs <- org.Hs.eg.db ### excel表格是从dbindel数据库中获取的:HCC1954 (indel.txt文件)文件<- read.csv("GSM721136.indel.txt"， sep="\t"， header = t) ids2 <- as.character(file$related_gene) #keytypes(hs) AnnotationDbi::mapIds(hs, keys = ids2, column='SYMBOL'， keytype='ENTREZID') ##成功地将ENTREZID转换为符号###检索所有ENTREZID的TSS(我们没有这些基因id的链信息)库(txdb . hsapiens . ucsc .hg19. knowngene)库(GenomicRanges) txdb <- txdb . hsapiens . ucsc .hg19. hg19. hg19. hg19。knownGene keytypes(txdb) [1] "CDSID" "CDSNAME" "EXONID" "EXONNAME" "GENEID" "TXID" "TXID" "TXID" "TXID" "TXID" "CDSNAME" "CDSSTART" "CDSSTRAND" "EXONCHROM" [8] "EXONEND" "EXONID" "EXONRANK" "EXONSTART" "EXONSTRAND" "GENEID" [15] "TXCHROM" "TXEND" "TXID" "TXNAME" "TXSTART" "TXSTRAND" "TXTYPE" select(txdb, keys=ids2, keytype = "GENEID"， columns=c("CDSCHROM"，"CDSSTART"，"CDSSTRAND"))

我想确认一下，这是否是获得EntrezIDs的TSS站点的正确方法。CDSstart可以作为TSS吗?如果是，为什么一个entrezID会有多个结果?mapID函数也可以用于此目的，因为它会给出以下错误。

如果方法是错误的，你能告诉我正确的包/方法/职位，因为我几乎迷失在这里?

mapIds_base(x，键，列，键类型，…， multiVals = multiVals): mapIds只能使用一个列。

我不认为CDS开始一般可以作为TSS，只有在成绩单没有5' UTR的情况下才可以。

以下是如何提取UCSC已知基因hg38轨道中所有转录本的TSS:

1.从UCSC knownGene轨道提取hg38的所有转录本:

library(TxDb.Hsapiens.UCSC.hg38. knowngene) tx <- transcripts(TxDb.Hsapiens.UCSC.hg38. txt)knownGene列= c(“tx_name”、“gene_id”)与247541年)tx #农庄对象范围和2元数据列:# seqnames范围链| tx_name gene_id # < Rle > < IRanges > < Rle > | <人物> < CharacterList > # [1] chr1 11869 - 14409 + | ENST00000456328.2 100287102 # [2] chr1 12010 - 13670 + | ENST00000450305.2 100287102 # [3] chr1 29554 - 31097 + | ENST00000473358.1 # [4] chr1 30267 - 31109 + | ENST00000469289.1 # [5] chr1 30366 - 30503 + | ENST00000607096.1 100302278  # ... ... ... ... . ... ...# [247537] chrUn_GL000220v1 155997-156149 + | enst00000620272v1 380608-380726 + | ENST00000620265.1 # [247539] chrUn_KI270442v1 217250-217401 - | ENST00000611690.1 # [247540] chrUn_KI270744v1 51009-51114 - | ENST00000616830.1 # [247541] chrUn_KI270750v1 148668-148843 + | ENST00000612925.1 # ------- # seqinfo: hg38基因组595序列(1个循环)

的gene_idmetadata列是一个CharacterList对象，将每个转录本链接到0或1个基因:

mcols(tx)$gene_id #字符长度247541 #[[1]]100287102 #[[2]]100287102 #[[3]]字符(0)#[[4]]字符(0)#[[5]]100302278 #[[6]]字符(0)#[[7]]字符(0)#[[8]]字符(0)#[[9]]字符(0)#[[10]]79501 #…# <247531更多元素>表(长度(mcols(tx)$gene_id)) # 01 # 50159 197382

实际上，将这个映射表示为字符向量可能会更方便:

mcols (tx)美元gene_id < as.character (mcols (tx) gene_id美元)tx #农庄组织与247541年对象范围和2元数据列:# seqnames范围链| tx_name gene_id # < Rle > < IRanges > < Rle > | <人物> <人物> # [1]chr1 11869 + | ENST00000456328.2 100287102 # [2] chr1 12010 + | ENST00000450305.2 100287102 # [3] chr1 29554 + | ENST00000473358.1 < NA > # [4] chr1 30267 + | ENST00000469289.1 < NA > # [5] chr1 30366 + | ENST00000607096.1 100302278  # ... ... ... ... . ... ...# [247537] chrUn_GL000220v1 155997 + | ENST00000619779.1 109864274 # [247538] chrUn_KI270442v1 380608 + | ENST00000620265.1  # [247539] chrUn_KI270442v1 217401 - | ENST00000611690.1  # [247540] chrUn_KI270744v1 51114 - | ENST00000616830.1  # [247541] chrUn_KI270750v1 148668 + | ENST00000612925.1  # ------- # seqinfo: hg38基因组的595个序列(1个循环)

2.将基因名称/符号添加到'tx'作为另一个元数据列:

TxDb.Hsapiens.UCSC.hg38中没有保存基因名称/符号。knownGene包，所以我们不能做成绩单(TxDb.Hsapiens.UCSC.hg38。knownGene, columns=c("tx_name"， "gene_id"， "SYMBOL"))来获得这一信息。相反，我们从org.Hs.eg.db包中获取它:

library(org. hs . exe .db) ENTREZID2SYMBOL <- select(org. hs . exe .db, mcols(tx)$gene_id, c("ENTREZID"， "SYMBOL")) # 'select()'返回many:1键和列之间的映射stopifnot(相同(ENTREZID2SYMBOL$ENTREZID, mcols(tx)$gene_id)) mcols(tx)$SYMBOL <- ENTREZID2SYMBOL$SYMBOL .

请注意，上面的另一种选择是使用Homo。sapiens包，它捆绑TxDb.Hsapiens.UCSC.hg19。knownGene和org.Hs.eg.db一起:

库(Homo.sapiens)。#包含GODb对象:GO.db ##包含关于:基因本体的数据##包含OrgDb对象:org.Hs.eg.db ##关于:智人的基因数据##分类Id: 9606 ##包含TxDb对象:TxDb. hsapiens . ucsc .hg19。knownGene ##关于智人的转录组数据##基于基因组:hg19 ## OrgDb基因id ENTREZID映射到TxDb基因id GENEID。成绩单(Homo。sapiens, columns=c("ENTREZID"， "SYMBOL")) # 'select()'返回键和列之间的1:1映射# GRanges对象有82960个范围和2个元数据列:# seqnames ranges strand b| SYMBOL ENTREZID #    |   # [1] chr1 11874-14409 + | DDX11L1 100287102 # [3] chr1 11874-14409 + | DDX11L1 100287102 # [4] chr1 69091-70008 + | OR4F5 79501 # [5] chr1 321084-321115 + |   # ... ... ... ... . ... ...# [82956] chrUn_gl000237 1-2686 - |   # [82957] chrUn_gl000241 20433-36875 - |   # [82958] chrUn_gl000243 11501-11530 + |   # [82959] chrUn_gl000243 13608-13637 + |   # [82960] chrUn_gl000247 5787-5816 - |   # ------- # seqinfo: hg19基因组93个序列(1个循环)

注意Homo的原因。sapiens的转录本比TxDb.Hsapiens.UCSC.hg38少。knownGene是默认的Homo。Sapiens基于hg19而不是hg38。可以通过以下方法进行更改:

TxDb(Homo.sapiens) <- TxDb. hsapiens . ucsc .hg38. knowngene

3.使用发起人()获得发起人或TSS:

我们使用上游= 0而且下游= 1以获得TSS(见?推动者在基因组特征包中获取更多信息):

tss <- promoters(tx, upstream=0, downstream=1) tss # GRanges对象，包含247541个范围和3个元数据列:# seqnames ranges strand | tx_name gene_id SYMBOL #    |    # [1] chr1 11869 + | ENST00000456328.2 100287102 DDX11L1 # [2] chr1 12010 + | ENST00000450305.2 100287102 DDX11L1 # [3] chr1 29554 + | ENST00000473358.1   # [4] chr1 30267 + | ENST00000469289.1   # [5] chr1 30366 + | ENST00000607096.1 100302278 MIR1302-2 # ... ... ... ... . ... ... ...# [247537] chrUn_GL000220v1 155997 + | ENST00000619779.1 109864274 RNA5-8SN4 # [247538] chrUn_KI270442v1 380608 + | ENST00000620265.1   # [247539] chrUn_KI270442v1 217401 - | ENST00000611690.1   # [247540] chrUn_KI270744v1 51114 - | ENST00000616830.1   # [247541] chrUn_KI270750v1 148668 + | ENST00000612925.1   # ------- # seqinfo:来自hg38基因组的595个序列(1个循环)

H。

TSS是转录本的第一个碱基，基因可以有多个转录本。NCBI基因ID属于基因，而不是转录本。因此，对于任何具有多个转录本的基因，你都应该期望有多个转录本起始位点。

不管怎样,使用选择和一个TxDb对象都可能是次优的。我会用Homo.sapiens包中。

> zz <- cdsBy(Homo。智人，"基因"，columns = c("ENTREZID"，"SYMBOL"))'select()'返回多个键和列之间的映射> zz GRangesList对象长度19809:$ ' 1 ' GRanges对象有8个范围和3个元数据列:seqnames范围链| cds_name符号< Rle > < IRanges > < Rle > | <人物> < CharacterList > [1] chr19 58858388 - 58858395 | < NA > A1BG [2] chr19 58858719 - 58859006 | < NA > A1BG [3] chr19 58861736 - 58862017 | < NA > A1BG [4] chr19 58862757 - 58863053 | < NA > A1BG [5] chr19 58863649 - 58863921 | < NA > A1BG [6] chr19 58864294 - 58864563 | < NA > A1BG [7] chr19 58864658 - 58864693 | < NA > A1BG [8] chr19 58864770 - 58864803 | < NA > A1BG ENTREZID < CharacterList > [1] 1 [2] [3] [4] 1 [5] 1 [6] 1 [7] 1 [8]1 ------- seqinfo: 93 sequences (1 circular) from hg19 genome $ ' 10 ' GRanges对象，1个范围和3个元数据列:seqnames ranges strand | cds_name SYMBOL    |   [1] chr8 18257514-18258386 + |  NAT2 ENTREZID  [1] 10 ------- seqinfo: 93 sequences (1 circular) from hg19 genome…< 19807 >元素# #调整得到TSS > unlist(调整(zz, 1))农庄对象与237874范围和3元数据列:seqnames范围链| cds_name象征ENTREZID < Rle > < IRanges > < Rle > | <人物> < CharacterList > < CharacterList > 1 chr19 58858395 - | < NA > A1BG 1 1 chr19 58859006 - | < NA > A1BG 1 1 chr19 58862017 - | < NA > A1BG 1 1 chr19 58863053 - | < NA > A1BG 1 1 chr19 58863921 - | < NA > A1BG 1  ... ... ... ... . ... ... ...9994 chr6 90575672 + |  CASP8AP2 9994 9994 chr6 90581008 + |  CASP8AP2 9994 9994 chr6 90581008 + |  CASP8AP2 9994 9994 chr6 90583522 + |  CASP8AP2 9994 9997 chr22 50962840 - |  SCO2 9997 ------- seqinfo: hg19基因组93条序列(1个循环)

真的非常谢谢詹姆斯·w·麦克唐纳而且Herve页面感谢你抽出时间回答我的问题。我想扩展这个线程，并添加如何在使用Ensembl注释时做到这一点。

library(biomaRt) library(regioneR) library(bsgenome . hapiens . ucsc .hg38) mart <- useDataset("hsapiens_gene_ensembl"， useMart("ensembl")) att <- listAttributes(mart) #以dataframe转录本的形式从biomaRt获取TSS信息<- getBM(attributes=c("chromosome_name"， "transcript_start"， "transcript_end"， "strand"， "ensembl_gene_id"，"gene_biotype"， "ensembl_transcript_id")， filters ="biotype"， values =c("protein_coding")，Mart = Mart) #将链的含义从1和-1更改为+和- transcripts$strand[transcripts$strand==-1] <- "-" transcripts$strand[transcripts$strand==1] <- "+" #将df转换为granges转录本。gr <- toGRanges(转录本)基因组(transcripts.gr) <-“hg38”#这很重要，因为启动子()函数使用链信息来正确选择负链链(transcripts.gr) <-转录本上基因的起始坐标。Gr $strand #从元数据列中删除strand属性，否则会抛出错误记录。Gr <- transcripts。gr[，-1] #可选:只保留规范染色体转录本:转录本。gr <- keepSeqlevels(transcripts. levels)gr, c(1:22，“X”，“Y”)，修剪。mode = "粗")#可选:将染色体转换为UCSC样式，如将1转换为chr1 seqlevelsStyle(transcripts.gr) <- "UCSC" #获取TSS启动子(transcripts.gr)。gr，上游=0，下游=1)#对于人类的启动子，通常TSS上游1500 bo和下游500 bp被认为是pr <-启动子(转录本。<- Read .csv("homo sapiens_genes. gr, upstream=1500, downstream=500) #阅读你想要获得的文件启动子序列的列表。tsv”,9 = " \ t”,头= t) #子集entrez基因的启动子对象id子< -子集(公关,公关ensembl_transcript_id % %文件single_exon_isoforms美元)#检索所有的启动子序列子集范围seq < -getSeq (BSgenome.Hsapiens.UCSC。hg38, sub) #使用entrezid名称更改序列名称(seq) <- sub$ensembl_transcript_id #以快速格式保存