《Modern Statistics for Modern Biology》Chapter 一： 离散数据模型的预测（1.1 – 1.3）
《Modern Statistics for Modern Biology》Chapter 一： 离散数据模型的预测（1.4 – 1.5）
《Modern Statistics for Modern Biology》Chapter 二: 统计建模（2.1-2.3）
《Modern Statistics for Modern Biology》Chapter 二: 统计建模（2.4-2.5）
《Modern Statistics for Modern Biology》Chapter 二 统计建模（2.6 – 2.7）
《Modern Statistics for Modern Biology》Chapter 二 统计建模（2.8 – 2.9）
从这章开始最开始记录一些markdown等小知识
$\hat{p}=\frac{1}{12}$：
掌握R语言中的apply函数族
卡方检验
Hardy-Weinberg equilibrium( 哈迪-温伯格平衡 )
带你理解beta分布
简单介绍一下R中的几种统计分布及常用模型
- 统计分布每一种分布有四个函数：
d――density（密度函数），p――分布函数，q――分位数函数，r――随机数函数。比如，正态分布的这四个函数为dnorm，pnorm，qnorm，rnorm。下面我们列出各分布后缀，前面加前缀d、p、q或r就构成函数名：norm：正态，t：t分布，f：F分布，chisq：卡方（包括非中心）unif：均匀，exp：指数，weibull：威布尔，gamma：伽玛，beta：贝塔lnorm：对数正态，logis：逻辑分布，cauchy：柯西，binom：二项分布，geom：几何分布，hyper：超几何，nbinom：负二项，pois：泊松signrank：符号秩，wilcox：秩和，tukey：学生化极差
如何预测一条序列是否含有CpG岛
2.10 示例：基因组中出现核苷酸模式
到目前为止，我们看到的例子集中在离散计数和分类数据的分布上。让我们来看一个距离分布的例子，它是准连续的。对基因组序列中特定基序实例之间距离分布的这一案例研究也将使我们能够探索生物导体中的特定基因组序列操作。
Biostrings包提供了对序列数据进行操作的工具，
DNAStringandDNAStringSet.是R中基础的数据结构。这些工具能让我对单条或者多条序列进行有效的操作。
library("Biostrings")
► 问题
通过浏览教程小片段，了解Biostring包中提供的一些有用的数据和函数。
| 函数 | 功能 |
|---|---|
| length | 返回序列长度 |
| names | 返回序列的名称 |
| [ | 从object中提取序列 |
| head; tail | 从object中提取第一条或者最后一条序列 |
| rev | 对序列进行反向排序 |
| width, nchar | 返回一个object中的所有序列的大小（比如碱基的数目） |
| ==, != | 两个object元素范围内的比较 |
| match, %in% | 匹配 |
| duplicated, unique | 去重复功能 |
| sort, order | 排序 |
| relist, split, extractList |
更多内容
Table 1:Low-level manipulation ofDNAStringSetandAAStringSetobjects.
Table 2:Counting / tabulating.
Table 3:Sequence transformation and editing.
Table 4:String matching / alignments.
► 解
- 第一行打印遗传密码信息，第二行返回
IUPAC核苷酸类型。第三行列出了Biostring软件包中所有可用的vignettes（有哪些参考说明书），第四行显示了一个特定的vignettes`（功能参数预览）。
> GENETIC_CODE
TTT TTC TTA TTG TCT TCC TCA TCG TAT TAC TAA TAG TGT TGC TGA TGG CTT CTC CTA CTG CCT CCC CCA CCG CAT CAC CAA CAG
"F" "F" "L" "L" "S" "S" "S" "S" "Y" "Y" "*" "*" "C" "C" "*" "W" "L" "L" "L" "L" "P" "P" "P" "P" "H" "H" "Q" "Q"
CGT CGC CGA CGG ATT ATC ATA ATG ACT ACC ACA ACG AAT AAC AAA AAG AGT AGC AGA AGG GTT GTC GTA GTG GCT GCC GCA GCG
"R" "R" "R" "R" "I" "I" "I" "M" "T" "T" "T" "T" "N" "N" "K" "K" "S" "S" "R" "R" "V" "V" "V" "V" "A" "A" "A" "A"
GAT GAC GAA GAG GGT GGC GGA GGG
"D" "D" "E" "E" "G" "G" "G" "G"
attr(,"alt_init_codons")
[1] "TTG" "CTG"
> IUPAC_CODE_MAP
A C G T M R W S Y K V H D B N
"A" "C" "G" "T" "AC" "AG" "AT" "CG" "CT" "GT" "ACG" "ACT" "AGT" "CGT" "ACGT"
> vignette(package = "Biostrings")
> vignette("BiostringsQuickOverview", package = "Biostrings")
-
BSgenome包提供了许多物种的基因组信息，您可以通过键入以下命令来访问包含整个基因组序列的数据包的名称
> library("BSgenome")
> ag = available.genomes()
> length(ag) ## 可以看到有91个基因组的信息
[1] 91
> ag[1:2]
[1] "BSgenome.Alyrata.JGI.v1" "BSgenome.Amellifera.BeeBase.assembly4"
- 接下来我们将探索
AGGAGGT在E.coli（大肠杆菌）中出现。我们使用一个特定菌株的基因组序列(Escherichia coli str. K12 substr.DH10B)，其NCBI登录号为NC_010473。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BSgenome.Ecoli.NCBI.20080805", version = "3.8")
library("BSgenome.Ecoli.NCBI.20080805")
> Ecoli
E. coli genome:
# organism: Escherichia coli (E. coli)
# provider: NCBI
# provider version: 2008/08/05
# release date: NA
# release name: NA
# 13 sequences:
# NC_008253 NC_008563 NC_010468 NC_004431 NC_009801 NC_009800 NC_002655 NC_002695 NC_010498 NC_007946 NC_010473 NC_000913
# AC_000091
# (use 'seqnames()' to see all the sequence names, use the '$' or '[[' operator to access a given sequence)
> shineDalgarno = "AGGAGGT"
> ecoli = Ecoli$NC_010473
- 我们可以使用
countPattern函数计算宽度为50000的窗口中出现的pattern。
> window = 50000
> starts = seq(1, length(ecoli) - window, by = window)
> ends = starts + window - 1
> numMatches = vapply(seq_along(starts), function(i) {
+ countPattern(shineDalgarno, ecoli[starts[i]:ends[i]],
+ max.mismatch = 0)
+ }, numeric(1))
> table(numMatches)
numMatches
0 1 2 3 4
48 32 8 3 2
> numMatches
[1] 0 1 0 1 4 0 0 0 1 0 0 1 0 0 1 0 0 0 0 3 0 0 4 0 1 0 2 2 1 0 1 0 1 1 1 0 0 0 0 1 3 1 1 1 1 0 0 0 0 1 0 0 1 0 2 1 2 1 1 0 0 0 0
[64] 0 1 0 1 1 0 1 2 0 0 0 0 0 0 0 0 1 1 0 3 1 0 1 2 1 2 1 2 0 1
> str(starts)
num [1:93] 1 50001 100001 150001 200001 ...
> str(numMatches)
num [1:93] 0 1 0 1 4 0 0 0 1 0 ...
►问题 What distribution might this table fit ?
► 解
- 我们可以使用
matchPattern函数检查AGGAGGT匹配到的位置信息。
> sdMatches = matchPattern(shineDalgarno, ecoli, max.mismatch = 0)
> sdMatches
Views on a 4686137-letter DNAString subject
subject: AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC...TGCTGCATGATATTGAAAAAAATATCACCAAATAAAAAACGCCTTAGTAAGTATTTTTC
views:
start end width
[1] 56593 56599 7 [AGGAGGT]
[2] 199644 199650 7 [AGGAGGT]
[3] 202176 202182 7 [AGGAGGT]
[4] 214433 214439 7 [AGGAGGT]
[5] 217429 217435 7 [AGGAGGT]
... ... ... ... ...
[61] 4438786 4438792 7 [AGGAGGT]
[62] 4498085 4498091 7 [AGGAGGT]
[63] 4536658 4536664 7 [AGGAGGT]
[64] 4546821 4546827 7 [AGGAGGT]
[65] 4611626 4611632 7 [AGGAGGT]
- 可以在R命令行中键入
sdMatches以获取此对象的summary。它包含所有65个模式匹配的位置，表示为一组所谓的原始序列视图。它们之间的距离是多少？
> betweenmotifs = gaps(sdMatches)
> betweenmotifs
Views on a 4686137-letter DNAString subject
subject: AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC...TGCTGCATGATATTGAAAAAAATATCACCAAATAAAAAACGCCTTAGTAAGTATTTTTC
views:
start end width
[1] 1 56592 56592 [AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAA...TGAATAATTTCGCTACCTACCACGCTCCAGGTGTCAGAACCCGGCAGAC]
[2] 56600 199643 143044 [AAAGCTACCGTTATCCAGAATGGTATAGCGGTTTTCGCCACGTTGTTTC...ACATGTTATGGGGCTATAGCTCAGCTGGGAGAGCGCCTGCTTTGCACGC]
[3] 199651 202175 2525 [CTGCGGTTCGATCCCGCATAGCTCCACCATCTCTGTAGTGGTTAAATAA...TAAAGAGTAACGGAGGAGCACGAAGGTTGGCTAATCCTGGTCGGACATC]
[4] 202183 214432 12250 [TAGTGCAATGGCATAAGCCAGCTTGACTGCGAGCGTGACGGCGCGAGCA...GATATTGGAAATATCTGATTTGCAAATTATCGTGTTATCGCCAGGCTTT]
[5] 214440 217428 2989 [TAATAACATGGGCAGGATAAGCTCGGGAGGAATGATGTTTAAGGCAATA...CCGTAGCGAGAATACTCAAAATCATCATAACGAAAAGCCCCTTACTTGT]
... ... ... ... ...
[62] 4438793 4498084 59292 [GGATTTAATCACGGTAACATTTTGCTGGCTTAAAGCATCTGCCAGACGC...GAGCAACAGGCGGCAGAGCAAGGTTGGGAGTCATTGCATCGTCAACTTC]
[63] 4498092 4536657 38566 [AGATCCGGTTGCGGCAGCAAGGATTCATCCAAATGATCCACAAAGGCTT...AGGATATAGTAAGCAAATTAAAACGGGGTAAATTTGAACTCCAAATACC]
[64] 4536665 4546820 10156 [GGAATTAAAGAATGCGATGGATGGTATATTTACGAAAAAATAATTGATG...GGGGTAAAACGTTTCCTGTAGCACCGTGAGTTATACTTTGTATAACTTA]
[65] 4546828 4611625 64798 [GCAGATGCGTATTACCATAAAAAGATGGGGGAACAGTGCAGGTATGGTC...ATGGGGCACCGACACACCGAGCGTGGTGGCGCGCGCATCGCCGAGTGCA]
[66] 4611633 4686137 74505 [CGAGATCGCGGCAAAAACTCAGGCTCAGCGGCAGAAATAAAATCATCAG...TATTGAAAAAAATATCACCAAATAAAAAACGCCTTAGTAAGTATTTTTC]
- 因此，这些实际上是66个互补的区域。现在，让我们找到一个模型，分布之间的差距大小的主题。如果模体发生在随机位置，我们期望间隙长度服从指数分布。下面的代码(其输出如
图2.23所示)评估了这一假设。如果指数分布是一个很好的拟合，点应该大致在一条直线上。指数分布只有一个参数，即速率，并给出了与数据估计相对应的斜率直线。
> library("Renext")
载入需要的程辑包：evd
Warning messages:
1: 程辑包‘Renext’是用R版本3.5.2 来建造的
2: 程辑包‘evd’是用R版本3.5.2 来建造的
> expplot(width(betweenmotifs), rate = 1/mean(width(betweenmotifs)),
+ labels = "fit")
图 2.23
2.10.1在依赖关系的情况下建模（Modeling in the case of dependencies）
正如我们在第2.8节中看到的，核苷酸序列通常是依赖的：
在给定的位置发现某个核苷酸的概率往往取决于周围的序列。在这里，我们将使用马尔可夫链进行依赖建模。我们将研究人类基因组第8号染色体的区域，并试图发现被称为CpG岛的区域与其他区域之间的差异。我们使用数据(由Irizarry，Wu，和Feinberg(2009)生成)，这些数据告诉我们
CpG的起点和终点在基因组中的什么位置，并观察核苷酸的频率和‘CG’、‘CT’、‘CA’、‘CC’的频率。因此，我们可以问，核苷酸的出现之间是否存在依赖关系，如果存在，如何对它们进行建模。
> library("BSgenome.Hsapiens.UCSC.hg19")
> chr8 = Hsapiens$chr8
> setwd("e://compute language/NGS/Modern Statistics for Modern Biology/data/")
> CpGtab = read.table("model-based-cpg-islands-hg19.txt",header = TRUE)
> nrow(CpGtab)
[1] 65699
> head(CpGtab)
chr start end length CpGcount GCcontent pctGC obsExp
1 chr10 93098 93818 721 32 403 0.559 0.572
2 chr10 94002 94165 164 12 97 0.591 0.841
3 chr10 94527 95302 776 65 538 0.693 0.702
4 chr10 119652 120193 542 53 369 0.681 0.866
5 chr10 122133 122621 489 51 339 0.693 0.880
6 chr10 180265 180720 456 32 256 0.561 0.893
> irCpG = with(dplyr::filter(CpGtab, chr == "chr8"), IRanges(start = start, end = end))
> head(irCpG)
IRanges object with 6 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 32437 33708 1272
[2] 40354 40656 303
[3] 44536 46203 1668
[4] 99457 100721 1265
[5] 155954 156761 808
[6] 179033 179756 724
- 在上面的行中，我们将数据集
CpGtab的子集(筛选)仅设为第8染色体，然后创建一个IRange对象，该对象的开始和结束位置由该数据集中命名相同的列定义。在IRange函数调用(它从其参数构造对象)中，第一个开始是函数的参数名称，第二个开始是作为filter的输出在数据集中获得的列；对于end也是类似的。IRange是一种通用的数学间隔容器。IRange中的“I”代表“间隔”，Granges中的“G”代表“基因组”。
> grCpG = GRanges(ranges = irCpG, seqnames = "chr8", strand = "+")
> grCpG
GRanges object with 2855 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr8 32437-33708 +
[2] chr8 40354-40656 +
[3] chr8 44536-46203 +
[4] chr8 99457-100721 +
[5] chr8 155954-156761 +
... ... ... ...
[2851] chr8 146126089-146128165 +
[2852] chr8 146175867-146176338 +
[2853] chr8 146228006-146228907 +
[2854] chr8 146232962-146233086 +
[2855] chr8 146276730-146278360 +
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
> genome(grCpG) = "hg19"
> library("Gviz")
载入需要的程辑包：grid
Warning message:
程辑包‘Gviz’是用R版本3.5.1 来建造的
> ideo = IdeogramTrack(genome = "hg19", chromosome = "chr8")
Warning messages:
1: In readChar(f, 1) : 在non-UTF-8 MBCS语言环境里只能读取字节
2: In readChar(f, 1) : 在non-UTF-8 MBCS语言环境里只能读取字节
> ideo
Ideogram track 'chr8' for chromosome 8 of the hg19 genome
> plotTracks(
+ list(GenomeAxisTrack(),
+ AnnotationTrack(grCpG, name = "CpG"), ideo),
+ from = 2200000, to = 5800000,
+ shape = "box", fill = "#006400", stacking = "dense")
- 我们现在定义所谓的染色体序列的观点，对应于
CpG岛，irCpG，和之间的区域(gap(IrCpG))。生成的对象CGIview和nonCGIview只包含坐标，而不包含序列本身(这些对象保存在对象Hsapiens$chr8中)，因此它们在存储方面是相当轻量级的。
> CGIview = Views(unmasked(Hsapiens$chr8), irCpG)
> CGIview
Views on a 146364022-letter DNAString subject
subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
views:
start end width
[1] 32437 33708 1272 [CGGGTGCCATCTCAGCAGCTCACGGTGTGGAAACTGCGACACTCACA...ATTGGAACTGTGCGTTTGCGGAAGAACACGGCGGGGGGGGCGGCGC]
[2] 40354 40656 303 [CGATTAGCGGGAGCGCAAAGCGAAGGGGCGGCCTGCGACGTTTTCCT...ACGAGGAGGGAGCCGCGCACGCGCGGTCGGCTCGGCGAGGAGCCGG]
[3] 44536 46203 1668 [CGGCGCTGCAGGAGAGGAGATGCCCAGGCCAGGCGGCCGGCGCACGT...TGTGATGTTGTGTTCTCGGTGTGAGTTCATGGGTGTGACGGGGCGT]
[4] 99457 100721 1265 [CGGCTGGCCGGGCAGGGGGGCTGACCCCCCCCACCTCCCTCCTGGAC...GCTCCGCTGCCCGCTGAACTCCATCCTCCCGGCGGTCGGGCGGCGG]
[5] 155954 156761 808 [CGGGAAAGATTTTATTCACCGTCGATGCGGCCCCGAGTTGTCCCAAA...GGCACGGCAGGGCTCTCTTGCTCGCAGTATACTGGCGGCACGCCGC]
... ... ... ... ...
[2851] 146126089 146128165 2077 [CGGCAAGTAGCACACTCGACCAACTGCTGAAACGAAAGCAGCTCGTT...AAAGAGAAGGCTTTATTATTTCAAGGTCGACGGCCAAGGAGACCGG]
[2852] 146175867 146176338 472 [CGACAGGCGGGAACGCCACCAGCCTGTGGGCGGAGTCCTGGCTCGGC...TATCTGGGGCAGCGTAGCCTCTCTAGATGGCGGGAGGCCGGGGCGA]
[2853] 146228006 146228907 902 [CGCCACCCTCGCCCTGCCCGCGCCTCGGCCAACACAAAGCGGGCGCT...GGGAGAGGCCAGGTGCGACCCGGGCAGGTGAGGAGCACGGGCCCGG]
[2854] 146232962 146233086 125 [CGCCTCCAGCTGCGCAAGGGCAGGTGCGGGCAAGGAAGTCCCACTGG...TGGGTGAGGGTAGAGAGGATGCCTCAGTGCGTCCGGACACCACCGA]
[2855] 146276730 146278360 1631 [CGTCTCAGGAAACACCGCTTTCCACTCCTGTGTCCCCGAAAGGAATC...ATCCACGGGTGCGAATCCACGGGCACGTGTGATCTGGGGTCCGCGG]
> NonCGIview = Views(unmasked(Hsapiens$chr8), gaps(irCpG))
> NonCGIview
Views on a 146364022-letter DNAString subject
subject: NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN...NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
views:
start end width
[1] 33709 40353 6645 [TCTTCTCATCACGTGCTCTCAGGGGCCACGATGTCAACATGCCTCA...AATGCTGGGTGCCCGCAAGAACGCGGAGCAGCAGGCACTCATTCC]
[2] 40657 44535 3879 [TCTCCAAGTGCCGCCAGCTGCGGGATTTCCTCTGCAAAAGACAAAC...CCAGGCCTGGCGGCCGGCGCACGCGGGTTCTCTGTGGCCAGCAGG]
[3] 46204 99456 53253 [GTGCTGTGTGAGAACATGTGTGTAGTGTTCACATGTCCTCTGTGCG...AGTAGGGGCGGCCGGGCAGAGGCGCCCCTCACCTCCCGGACGGGG]
[4] 100722 155953 55232 [CGGCTGCCCTAAGTGATCTTTTTAAGTAAAGGAGAAACTCACCAGG...AACCGTTCTAACTGGTCTCTGACCTTGATTATTAACGGCTGCAAC]
[5] 156762 179032 22271 [CTGCTGGCAGCTAGGGACATTGCAGGGCCCTCTTGCTCACATTGTA...CCATGTCTGCATCCCGGGCGGTATGTACGCTCTGAGAGATACATG]
... ... ... ... ...
[2850] 146125886 146126088 203 [GGGGGATGACTTTTATATATTTCTGACATGCCTGCATTAGAAGAAG...AGAAAAGACAGCATCGACAGCAAAATGGGGAGAGGATAAGACACC]
[2851] 146128166 146175866 47701 [AGGCATGGCTCAAATGTGGCTCCCACATGTGGGCTTAGGTCAAAGT...CCACCCCTTCACCCCGTCGCTCCGCCCGAGCTCGAATGGTGGCAC]
[2852] 146176339 146228005 51667 [CTCTGTGGTTCAGCCCCAGACCTGGTCCCACTAGCCCTGAGGCCAG...ACAGCCCGTAGACCTGCACGCCCCGCTGCCTCCGCTCGTGCCTCA]
[2853] 146228908 146232961 4054 [AGGGAGGACTTACGAGAATGGAAAGGAAAGCAGCCCGTCTTCAGTG...CTGAAGACAGCGCAGCTATCCGTAAGGATCACATGCGGCTTCCCT]
[2854] 146233087 146276729 43643 [GACCATAGTTCTCCAAGGGCACCTAGAGTGGATCCAGCTGAGTTCA...CCTTCGAGACCTTGCAGACCAGGCAGCTAACACTTCCCACCCCAC]
- 我们利用这些数据计算了
CpG岛和非CpG岛的跃迁计数。
> seqCGI = as(CGIview, "DNAStringSet")
> seqCGI
A DNAStringSet instance of length 2855
width seq
[1] 1272 CGGGTGCCATCTCAGCAGCTCACGGTGTGGAAACTGCGACACTCACACGGGTGCCATC...TCTCTGATTACATTGGAACTGTGCGTTTGCGGAAGAACACGGCGGGGGGGGCGGCGC
[2] 303 CGATTAGCGGGAGCGCAAAGCGAAGGGGCGGCCTGCGACGTTTTCCTGTAAAGTTGGG...GCTCAGCACGGACGAGGAGGGAGCCGCGCACGCGCGGTCGGCTCGGCGAGGAGCCGG
[3] 1668 CGGCGCTGCAGGAGAGGAGATGCCCAGGCCAGGCGGCCGGCGCACGTGGGTTCTCTGT...GAGTCCCTGTGTGTGATGTTGTGTTCTCGGTGTGAGTTCATGGGTGTGACGGGGCGT
[4] 1265 CGGCTGGCCGGGCAGGGGGGCTGACCCCCCCCACCTCCCTCCTGGACGGGGCGGCTGG...AGACTGAGACAGCTCCGCTGCCCGCTGAACTCCATCCTCCCGGCGGTCGGGCGGCGG
[5] 808 CGGGAAAGATTTTATTCACCGTCGATGCGGCCCCGAGTTGTCCCAAAGCCAGGCAGTG...GCTGGCGACGGGGCACGGCAGGGCTCTCTTGCTCGCAGTATACTGGCGGCACGCCGC
... ... ...
[2851] 2077 CGGCAAGTAGCACACTCGACCAACTGCTGAAACGAAAGCAGCTCGTTTTCGTCTTGAC...AGTTTTTGAGAAAAGAGAAGGCTTTATTATTTCAAGGTCGACGGCCAAGGAGACCGG
[2852] 472 CGACAGGCGGGAACGCCACCAGCCTGTGGGCGGAGTCCTGGCTCGGCCCCTCGCGGCG...GCCGCTCGTTCTATCTGGGGCAGCGTAGCCTCTCTAGATGGCGGGAGGCCGGGGCGA
[2853] 902 CGCCACCCTCGCCCTGCCCGCGCCTCGGCCAACACAAAGCGGGCGCTCAGCACCCCGC...AGCGGAGTGAGGGGAGAGGCCAGGTGCGACCCGGGCAGGTGAGGAGCACGGGCCCGG
[2854] 125 CGCCTCCAGCTGCGCAAGGGCAGGTGCGGGCAAGGAAGTCCCACTGGGGCTTTTGGGG...GCTCGGGCGAGTGGGTGAGGGTAGAGAGGATGCCTCAGTGCGTCCGGACACCACCGA
[2855] 1631 CGTCTCAGGAAACACCGCTTTCCACTCCTGTGTCCCCGAAAGGAATCATCAGACATCT...GGAAGGGTGGAATCCACGGGTGCGAATCCACGGGCACGTGTGATCTGGGGTCCGCGG
> seqNonCGI = as(NonCGIview, "DNAStringSet")
> dinucCpG = sapply(seqCGI, dinucleotideFrequency)
> dinucNonCpG = sapply(seqNonCGI, dinucleotideFrequency)
> dinucNonCpG[, 1]
AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT
389 351 400 436 498 560 112 603 359 336 403 336 330 527 519 485
> dinucCpG[, 1]
AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT
49 127 55 8 132 44 116 112 54 129 107 100 4 104 112 18
> NonICounts = rowSums(dinucNonCpG)
> IslCounts = rowSums(dinucCpG)
> NonICounts
AA AC AG AT CA CC CG CT GA GC GG GT TA TC
14223322 7129981 9795572 11291315 10241505 6931732 1174927 9779692 8372468 5706958 6934755 7112531 9602902 8359398
TG TT
10221420 14197986
> IslCounts
AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT
64103 83109 122131 44000 113346 193847 133549 122377 105186 177326 194517 86752 31210 106738 114592 65861
- For a four state Markov chain as we have, we define the transition matrix as a matrix where the rows are the from state and the columns are the to state. (这里放原文比较好理解)
> TI = matrix( IslCounts, ncol = 4, byrow = TRUE)
> TI
[,1] [,2] [,3] [,4]
[1,] 64103 83109 122131 44000
[2,] 113346 193847 133549 122377
[3,] 105186 177326 194517 86752
[4,] 31210 106738 114592 65861
> TnI = matrix(NonICounts, ncol = 4, byrow = TRUE)
> dimnames(TI) = dimnames(TnI) =
+ list(c("A", "C", "G", "T"), c("A", "C", "G", "T"))
> TI
A C G T
A 64103 83109 122131 44000
C 113346 193847 133549 122377
G 105186 177326 194517 86752
T 31210 106738 114592 65861
- 我们使用每种类型的跃迁次数的计数来计算频率，并将它们放入两个矩阵中 。
> TI
A C G T
A 64103 83109 122131 44000
C 113346 193847 133549 122377
G 105186 177326 194517 86752
T 31210 106738 114592 65861
> MI = TI /rowSums(TI)
> MI
A C G T
A 0.20457773 0.2652333 0.3897678 0.1404212
C 0.20128250 0.3442381 0.2371595 0.2173200
G 0.18657245 0.3145299 0.3450223 0.1538754
T 0.09802105 0.3352314 0.3598984 0.2068492
> MN = TnI / rowSums(TnI)
> MN
A C G T
A 0.3351380 0.1680007 0.23080886 0.2660524
C 0.3641054 0.2464366 0.04177094 0.3476871
G 0.2976696 0.2029017 0.24655406 0.2528746
T 0.2265813 0.1972407 0.24117528 0.3350027
► 问题
- 这意味着不同行中的转换是否不同？例如，
► 解
- 转变是不同的。例如，对于CpG岛(MI)转换矩阵，从C到A的转换和从T到A的转换看起来非常不同(0.201和0.098)。
► 问题
- CpG岛上不同核苷酸的相对频率是否不同？
► 解
> freqIsl = alphabetFrequency(seqCGI, baseOnly = TRUE, collapse = TRUE)[1:4]
> freqIsl / sum(freqIsl)
A C G T
0.1781693 0.3201109 0.3206298 0.1810901
> freqNon = alphabetFrequency(seqNonCGI, baseOnly = TRUE, collapse = TRUE)[1:4]
> freqNon / sum(freqNon)
A C G T
0.3008292 0.1993832 0.1993737 0.3004139
- 从结果可以看出，在
CpG岛中，C和G的频率大概是0.32，在非-CpG岛区，有A和T的频率为0.30。
► 问题
- 我们如何利用这些差异来决定一个给定的序列是否来自CpG岛？
► 解
使用
χ2-平方（卡方检验）统计量来比较观察到的频率与freqIsI和freqNon频率之间的差异。对于较短的序列，这可能不够敏感，下面给出一个更敏感的方法。给定一个序列，我们不知道它是否在一个CpG岛上，我们可以问它与其他地方相比，它属于一个CpG岛的概率有多大。我们根据所谓的优势比来计算分数。让我们举个例子：假设我们的序列x是
ACGTTATACTACG，我们想要决定它是否来自CpG岛。假设序列来自CpG岛，如果我们将该序列建模为一阶
Markov链，则可以这样写：
我们将把这个概率与
非-CpG岛的概率进行比较。正如我们在上面看到的，这些概率往往是非常不同的。我们将取它们的比值，看看它是大于还是小于1。这些主权将是许多小term的产物，并将变得非常小。我们可以通过取对数来解决这个问题。
这是对数似然比分数。为了加快计算速度，我们计算了对数比
,然后对相关的数据进行汇总得到我们的分数。
> alpha = log((freqIsl/sum(freqIsl)) / (freqNon/sum(freqNon)))
> alpha
A C G T
-0.5238084 0.4734387 0.4751059 -0.5061665
> beta = log(MI / MN)
> beta
A C G T
A -0.4935945 0.4566418 0.5239611 -0.6390469
C -0.5927341 0.3342289 1.7365319 -0.4699321
G -0.4671646 0.4383576 0.3360277 -0.4967508
T -0.8379216 0.5303960 0.4002977 -0.4821485
> x = "ACGTTATACTACG"
> scorefun = function(x) {
+ s = unlist(strsplit(x, ""))
+ score = alpha[s[1]]
+ if (length(s) >= 2)
+ for (j in 2:length(s))
+ score = score + beta[s[j-1], s[j]]
+ score
+ }
> scorefun(x)
A
-0.2824623
- 在下面的代码中，我们从
seqCGI对象中的2855个序列中选择长度为len=100的序列，然后从seqNonCGI对象中的2854个序列中选择(每个序列都是一个DNAStringSet)。在generateRandomScore函数的前三行中，我们删除了包含除A、C、T、G以外的任何字母的序列，例如"."(用于未定义核苷酸的字符)。在其余的序列中，我们以与其长度减去len成正比的概率进行采样，然后从中提取长度len的序列。对子序列的起始点进行均匀采样，同时满足子序列必须满足的约束条件。
> generateRandomScores = function(s, len = 100, B = 1000) {
+ alphFreq = alphabetFrequency(s)
+ isGoodSeq = rowSums(alphFreq[, 5:ncol(alphFreq)]) == 0
+ s = s[isGoodSeq]
+ slen = sapply(s, length)
+ prob = pmax(slen - len, 0)
+ prob = prob / sum(prob)
+ idx = sample(length(s), B, replace = TRUE, prob = prob)
+ ssmp = s[idx]
+ start = sapply(ssmp, function(x) sample(length(x) - len, 1))
+ scores = sapply(seq_len(B), function(i)
+ scorefun(as.character(ssmp[[i]][start[i]+(1:len)]))
+ )
+ scores / len
+ }
> scoresCGI = generateRandomScores(seqCGI)
> scoresNonCGI = generateRandomScores(seqNonCGI)
> br = seq(-0.6, 0.7, length.out = 50)
> h1 = hist(scoresCGI, breaks = br, plot = FALSE)
> h2 = hist(scoresNonCGI, breaks = br, plot = FALSE)
> plot(h1, col = rgb(0, 0, 1, 1/4), xlim = c(-0.5, 0.5), ylim=c(0,120))
> plot(h2, col = rgb(1, 0, 0, 1/4), add = TRUE)
图 2.25
2.11 Summary of this chapter
Goodness of fit 我们使用了不同的可视化，并展示了如何运行模拟实验来检验我们的数据是否可以用一个公平的四盒多项式模型进行拟合。我们遇到了
卡方检验，看到了如何使用QQ图将模拟和理论进行比较。Estimation 我们解释了
最大似然和贝叶斯估计过程。举例说明了这些方法，包括核苷酸模式发现和单倍型估计。Prior and posterior distributions 在评估以前研究过的类型的数据(如单倍型)时，计算数据的后验分布可能是有益的。这使我们能够通过简单的计算，在决策过程中加入不确定性。只要有足够的数据，先验的选择对结果影响不大。
CpG islands and Markov chains 我们了解了如何通过
马尔可夫链转移来模拟DNA序列中的依赖关系。我们利用这一点来建立基于似然比的分数，这样我们就可以看到长DNA序列是否来自CpG岛。当我们制作分数直方图时，我们在图 2.25中看到了一个值得注意的特征：它似乎是由两个部分组成的。这种双峰是我们第一次遇到混合物，它们是第四章的主题。这是在一些训练数据上建立模型的第一个例子：我们知道序列在CpG岛上，我们以后可以用这些序列对新的数据进行分类。我们将在
第12章中制定一种更完整的方法来实现这一点。
References
Cleveland, William S. 1988. The Collected Works of John W. Tukey: Graphics 1965-1985. Vol. 5. CRC Press.
Rice, John. 2006. Mathematical Statistics and Data Analysis. Cengage Learning.
Elson, D, and E Chargaff. 1952. “On the Desoxyribonucleic Acid Content of Sea Urchin Gametes.” Experientia 8 (4). Springer: 143–45.
Mourant, AE, Ada Kopec, and K Domaniewska-Sobczak. 1976. “The Distribution of the Human Blood Groups 2nd Edition.” Oxford University Press London.
Finetti, Bruno de. 1926. “Considerazioni Matematiche Sull’ereditarieta Mendeliana.” Metron 6: 3–41.
Cannings, Chris, and Anthony WF Edwards. 1968. “Natural Selection and the de Finetti Diagram.” Annals of Human Genetics 31 (4). Wiley Online Library: 421–28.
Love, Michael I, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2.” Gnome Biology 15 (12): 1–21.
Irizarry, Rafael A, Hao Wu, and Andrew P Feinberg. 2009. “A Species-Generalized Probabilistic Model-Based Definition of Cpg Islands.” Mammalian Genome 20 (9-10). Springer: 674–80.
Durbin, Richard, Sean Eddy, Anders Krogh, and Graeme Mitchison. 1998. Biological Sequence Analysis. Cambridge University Press.
Freedman, David, Robert Pisani, and Roger Purves. 1997. Statistics. New York, NY: WW Norton.
Agresti, Alan. 2007. An Introduction to Categorical Data Analysis. John Wiley.
Grantham, Richard, Christian Gautier, Manolo Gouy, M Jacobzone, and R Mercier. 1981. “Codon Catalog Usage Is a Genome Strategy Modulated for Gene Expressivity.” Nucleic Acids Research 9 (1). Oxford Univ Press: 213–13.
Perrière, Guy, and Jean Thioulouse. 2002. “Use and Misuse of Correspondence Analysis in Codon Usage Studies.” Nucleic Acids Research 30 (20). Oxford Univ Press: 4548–55.
Robert, Christian, and George Casella. 2009. Introducing Monte Carlo Methods with R. Springer Science & Business Media.
Marin, Jean-Michel, and Christian Robert. 2007. Bayesian Core: A Practical Approach to Computational Bayesian Statistics. Springer Science & Business Media.
McElreath, Richard. 2015. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. Chapman; Hall/CRC.
