Bioconductor没想象的那么简单-part5

豆豆写于19.4.11-12
keywords：豆花380  2  🐑城
Bioconductor没想象的那么简单(part1)
Bioconductor没想象的那么简单(part2)
Bioconductor没想象的那么简单(part3)
Bioconductor没想象的那么简单-part4

AnnotationHub

主要存储了一些公共数据和注释信息

library(AnnotationHub)
options()$BioC_mirror ##查看使用bioconductor的默认镜像
# 检查是否存在新版本，先赋一个值，有更新则自动更新
ahub <- AnnotationHub() #目前最新版本是2018-10-24
# 看看hub中都包括什么内容
> names(mcols(ahub))# 看看hub中都包括什么内容
 [1] "title"              "dataprovider"       "species"           
 [4] "taxonomyid"         "genome"             "description"       
 [7] "coordinate_1_based" "maintainer"         "rdatadateadded"    
[10] "preparerclass"      "tags"               "rdataclass"        
[13] "rdatapath"          "sourceurl"          "sourcetype"   

公共数据提供者dataprovider

unique(ahub$dataprovider) #公共数据提供者
table(ahub$dataprovider)#查看每个提供者都有多少数据

其中包含了多少物种基因组？

> head(unique(ahub$species))
[1] "Homo sapiens"         "Vicugna pacos"        "Dasypus novemcinctus"
[4] "Otolemur garnettii"   "Papio hamadryas"      "Papio anubis"  
> length(unique(ahub$species)) #包含了多少物种基因组
[1] 2042

hub提供哪些类的信息呢？

> table(ahub$rdataclass) #提供哪些类的信息

                      AAStringSet                        BigWigFile 
                                1                             10247 
                           biopax                         ChainFile 
                                9                              1113 
                       data.frame data.frame, DNAStringSet, GRanges 
                               40                                 3 
                            EnsDb                           GRanges 
                              678                             20475 
                           igraph                     Inparanoid8Db 
                                1                               268 
                             list                            MSnSet 
                               18                                 1 
                         mzRident                           mzRpwiz 
                                1                                 1 
                            OrgDb                               Rle 
                             2002                              1852 
                 SQLiteConnection                        TwoBitFile 
                                1                              5558 
                             TxDb                           VcfFile 
                               75                                 8 

• 其中FaFile 就是fasta文件

• UCSC的bigwig格式存储了基因组信号值coverage情况，用每个坐标区间的测序深度来表示。它是wig格式的压缩后的带有索引的二进制格式[利用wigToBigWig工具]，速度更快。

关于测序深度，可能我们对BAM文件比较熟悉，利用samtools depth input.bam可以计算测序深度（基因组上每bp碱基覆盖到的reads数目），一般得到的结果文件是：第一列是染色体，第二列是比对到染色体上不同碱基的位置，第三列为覆盖第二列位点的reads数目。

既然如此，为什么还要wig/bigwig或bedgraph文件呢？

因为人类基因组3G，BAM文件是非常大的，这种往往造成后续的可视化分析不好进行，比如一个3G左右的BAM文件，其包含的coverage信息在bigwig中只有十几M，一般ChIP-seq的覆盖度、测序深度分布在IGV的可视化中，使用bigwig比较频繁。这里只建立基础的认知，具体数据格式及可视化日后再来。

• chainfile 的作用就是进行基因组不同版本之间的坐标转换，如hg17/hg18/hg19/hg38

• OrgDb是做富集分析时主要用到的

BioStrings & BSgenome

这个包主要是用来处理基因序列的

准备基因组可以利用AnnotationHub或者BSgenome，其实AnnotationHub的基因组资源更多，但这里为了全面性考虑，主要介绍BSgenome的方法

BSgenome可以提供一个标准的注释信息，主要针对模式生物
其中关于人的，包括了hg17、18、19、38
小tip：

• 为什么看到下载的基因组中有N，还基本连续排列？表示没有测序出来的碱基，可以是ATCG的任意碱基

• 然后标准的CDS序列是ATCG四个大写字母

• 为什么还有小写的碱基？表示虽然测序成功但是序列复杂度比较低，一般是内含子（在成熟mRNA形成之前会被切除掉）或UTR（非编码）区域或者重复区域 =》 soft mask序列

• 最后，如果基因结构有不明白的，推荐看 下 https://vip.biotrainee.com/d/112–/4

BiocManager::install("BSgenome", version = "3.8")
library(BSgenome)
available.genomes() #总共91个,现在下载hg38
BiocManager::install("BSgenome.Hsapiens.UCSC.hg38", version = "3.8")
library("BSgenome.Hsapiens.UCSC.hg38")
# ?BSgenome.Hsapiens.UCSC.hg38
hg38 <- BSgenome.Hsapiens.UCSC.hg38

BSgenome的hg38基因组中都有什么信息？

可以看到物种、提供者、版本、发布日期、发布标准名称、染色体名称

染色体除了熟知的22对常染色体+X+Y+线粒体的，还有像chrUn_KI…这种，表示可以测出来确定是人的序列但拼不回去，还有像这种chr4_GL000008v2_random表示确定是人的4号染色体的序列但拼不回去。想看全部的名称使用seqnames(hg38)

> hg38
Human genome:
# organism: Homo sapiens (Human)
# provider: UCSC
# provider version: hg38
# release date: Dec. 2013
# release name: Genome Reference Consortium GRCh38
# 455 sequences: 
#   chr1                    chr2                    chr3                   
#   chr4                    chr5                    chr6                   
#   chr7                    chr8                    chr9                   
#   chr10                   chr11                   chr12                  
#   chr13                   chr14                   chr15                  
#   ...                     ...                     ...                    
#   chrUn_KI270744v1        chrUn_KI270745v1        chrUn_KI270746v1       
#   chrUn_KI270747v1        chrUn_KI270748v1        chrUn_KI270749v1       
#   chrUn_KI270750v1        chrUn_KI270751v1        chrUn_KI270752v1       
#   chrUn_KI270753v1        chrUn_KI270754v1        chrUn_KI270755v1       
#   chrUn_KI270756v1        chrUn_KI270757v1                               
# (use 'seqnames()' to see all the sequence names, use the '$' or '[['
# operator to access a given sequence)

截取序列

seqlengths(hg38) # 获取上面seqnames各条染色体的长度
hg38$chr1 #或者 getSeq(hg38, "chr1")
# 可以像这样创建
str1 = getSeq(hg38, "chr1", start = 10001, end = 10030)
# 或者直接创建
# str1 = DNAString("TACCCTAACCCTAACTAACCCTAA")

小操作：使用ggplot2对各条染色体长度进行可视化

# 使用ggplot2对各条染色体长度进行可视化
if(F){
  hg38.len <- data.frame(chr_name = seqnames(hg38),
                         chr_len = as.numeric(seqlengths(hg38)))
  hg38.len <- hg38.len[1:25,]
  #为了避免按ASCII编码chr1、chr10、chr11...画图，我们先进行因子转化，保证画出来的是chr1、chr2、chr3...
  hg38.len$chr_name <- factor( hg38.len$chr_name,
                               levels = as.character( hg38.len$chr_name))
  library(ggplot2)
  ggplot(hg38.len, aes(x = chr_name, y = chr_len)) + 
    geom_bar(stat = "identity")

}

可以非常直观看到：人类最小的染色体是21号

然后进行简单的探索

> reverse(str1) #反向序列
  24-letter "DNAString" instance
seq: AATCCCAATCAATCCCAATCCCAT
> reverseComplement(str1) #反向互补序列
  24-letter "DNAString" instance
seq: TTAGGGTTAGTTAGGGTTAGGGTA
> str1[1:15] #取其中一部分
  15-letter "DNAString" instance
seq: TACCCTAACCCTAAC
> length(str1) #统计长度
[1] 24
> translate(str1) #CDS序列按照标准的密码子表翻译成氨基酸
  8-letter "AAString" instance
seq: YPNPN*P* #其中*表示stop codon
> alphabetFrequency(str1) #统计str1中每个碱基出现的次数
 A  C  G  T  M  R  W  S  Y  K  V  H  D  B  N  -  +  . 
 9 10  0  5  0  0  0  0  0  0  0  0  0  0  0  0  0  0 

可以看到上面统计str1中每个碱基出现的次数时除了常规的ATCGN，还有其他的一些字母，它们其实表示了简并碱基，或者可以看看国际标准的兼并碱基

> IUPAC_CODE_MAP
     A      C      G      T      M      R      W      S      Y      K      V 
   "A"    "C"    "G"    "T"   "AC"   "AG"   "AT"   "CG"   "CT"   "GT"  "ACG" 
     H      D      B      N 
 "ACT"  "AGT"  "CGT" "ACGT" 

统计GC信息：

> letterFrequency(str1,"GC") #统计GC信息
G|C 
 10 
> letterFrequency(str1,"GC",as.prob = T) #统计GC信息(百分比)
      G|C 
0.4166667 
> dinucleotideFrequency(str1,as.prob = T) #每个两两碱基组合百分比
       AA        AC        AG        AT        CA        CC        CG 
0.1739130 0.1739130 0.0000000 0.0000000 0.0000000 0.2608696 0.0000000 
       CT        GA        GC        GG        GT        TA        TC 
0.1739130 0.0000000 0.0000000 0.0000000 0.0000000 0.2173913 0.0000000 
       TG        TT 
0.0000000 0.0000000 
# trinucleotideFrequency(str1,as.prob = T) # 当然统计三个碱基组合也可以

取基因组指定染色体的一段序列

> str2 <- hg38[["chr1"]][10000:20000] #取基因组指定染色体的一段
> str2
  10001-letter "DNAString" instance
seq: NTAACCCTAACCCTAACCCTAACCCTAACCCTAAC...AACTGAGACTGGGGAGGGACAAAGGCTGCTCTGT
# 然后可以像上面一样统计GC含量，比如统计chr2的GC含量
> letterFrequency(hg38[["chr2"]],"GC",as.prob = T)
      G|C 
0.3995527