1介绍

病毒准种目前被定义为一个密切相关的病毒基因组的集合，它经历了一个持续的遗传变异过程，产生的变异之间的竞争，并在给定的环境中选择最适合的基因组。(Domingo, Sheldon和Perales 2012)

产生这种准种的高复制错误率是由于缺乏遗传校对机制，据估计，对于典型的高复制率的病毒，每一个可能的点突变和许多双突变都是在病毒的每一个复制周期中产生的，这些可能随时出现在种群中。

准种的复杂性可以解释或预测病毒的行为;因此，它具有明显的临床意义。我们通常感兴趣的是比较来自单个患者的序列样本之间或来自不同患者组的样本之间的病毒多样性指数。这些比较可以为患者的临床进展或特定治疗的适当性提供信息。(Gregori et al. 2016)(Gregori et al. 2014)

QSUtils是一个用于NGS获得的准种扩增子数据的包，但它也可以用于分析16S/18S核糖体为基础的宏基因组学或通过扩增子的肿瘤遗传多样性。

在本教程中，我们演示了包中提供的函数的使用，以探索和操作序列对齐、将读取转换为单倍型和频率、修复读取、交叉链单倍型和可视化单倍型对齐。

2安装包

BiocManager::install("QSutils")

##警告:当版本与当前版本相同或大于当前版本时，未安装包;使用## ' force = TRUE '重新安装:'QSutils'

库(QSutils)

3.数据

该包包含处理准种数据的函数，准种数据由单倍型及其频率的对齐定义。数据从fasta格式的文件中加载，其中每个序列的头描述了单倍型的ID及其在准种种群中的相应频率。这两条信息被一条竖线“|”隔开。当频率信息缺失时，每个序列都被认为是一次读取。

filepath <执行(“extdata”、“ToyData_10_50_1000。fna"， package="QSutils") cat(readLines(filepath)， sep =" \n")

# # > Hpl_0_0001 | 464 | 46.4 # # CACCCTGTGACCAGTGTGTGCGAGTTTGGTGCCCCGTTGGCATTGACTAT # # > Hpl_1_0001 | 62 | 6.2 # # CACTCTGTGACCAGTGTGTGCGAGTTTGGTGCCCCGTTGGCATTGACTAT # # > Hpl_1_0002 | 39 | 3.9 # # CACCCTGTGACCAGCGTGTGCGAGTTTGGTGCCCCGTTGGCATTGACTAT # # > Hpl_1_0003 | 2.7 | 27 # # CACCCTGTGACCAGTGTGTGCGAGTTTGGTGCCCCGTTGGCATTGACTAC # # > Hpl_2_0001 | 3.7 | 37 # # CACTCTGTGACCAGTGTGTGCGAGTTTGGTGCCCCGTTAGCATTGACTAT # # > Hpl_4_0001 | 16 | 1.6 # # CACTCGGTGACCAGTGTGCGTGAGTTTGGTGCCCCGTTGGCATTGACTAT # # > Hpl_5_0001 | 33 # # | 3.3## cactctgtgaccagtgtgcgctttagtgccctgccggcatcgactat ## >Hpl_9_0001|248|24.8 ## cactctgtgtgatcagtgtgcgctttagtgccctgccggcatcgactac ## >Hpl_10_0001|20|2 ## cactctgtgtgtgcagtgtgcgctttagtgccctgccggcaccactac

的ReadAmplSeqs函数从fasta文件加载数据并返回一个包含两个元素的列表:DNAStringSet工厂单倍型序列和计数向量，nr．

filepath <执行(“extdata”、“ToyData_10_50_1000。fna"， package="QSutils") lst <- ReadAmplSeqs(filepath,type="DNA"

## $nr ## [1] 464 62 39 27 37 16 33 54 248 20 ## ## $hseqs ## DNAStringSet长度为10的对象:## width seq names ## [1] 50 caccctgaccagtgtgtgcga…gtgccccgttggcatgactat Hpl_0_0001|464|46.4 ## [2] 50 CACTCTGTGACCAGTGTGTGCGA…GTGCCCCGTTGGCATTGACTAT Hpl_1_0001|62|6.2 ## [3] 50 CACCCTGTGACCAGCGTGTGCGA…GTGCCCCGTTGGCATTGACTAT Hpl_1_0002|39|3.9 ## [4] 50 CACCCTGTGACCAGTGTGTGCGA…gtgccccgttggcatgactac Hpl_1_0003|27|2.7 ## [5] 50 CACTCTGTGACCAGTGTGTGCGA…GTGCCCCGTTAGCATTGACTAT Hpl_2_0001|37|3.7 ## [6] 50 CACTCGGTGACCAGTGTGCGTGA…GTGCCCCGTTGGCATTGACTAT Hpl_4_0001|16|1.6 ## [7] 50 CACTCTGTGACCAGTGTGCGCGA…GTGCCCTGTCGGCATTGACTAT Hpl_5_0001|33|3.3 ## [8] 50 CACTCTGTGATCAGTGTGCGCGA…GTGCCCTGCCGGCATCGACTAT Hpl_8_0001|54|5.4 ## [9] 50 CACTCTGTGATCAGTGTGCGCGA…GTGCCCTGCCGGCATCGACTAC Hpl_9_0001|248|24.8 ## [10] 50 CACTCTGTGATCAGTGTGCGCGA…GTGCCCTGCCGGCACCGACTAC Hpl_10_0001 20 | 2 |

的GetQSData函数从fasta文件加载数据，删除相对丰度低于给定阈值的单倍型，并对剩余的单倍型进行排序，首先通过增加相对于主要单倍型的突变数量，然后通过减少突变数量内的频率。

filepath <执行(“extdata”、“ToyData_10_50_1000。fna"， package="QSutils") lstG <- GetQSData(filepath,min。pct= 2,type="DNA"

长度为9的DNAStringSet对象:## width seq names ## [1] 50 CACCCTGTGACCAGTGTGTGCGA…仙人掌tgtgaccagtgtgtgcga…ggtgccccgttggcatgactat Hpl_1_0001 ## [3] 50 caacctgtgaccagcgtgtgcga…ggtgccccgttggcatgactat Hpl_1_0002 ## [4] 50 CACCCTGTGACCAGTGTGTGCGA…ggtgccccgttggcatgactac Hpl_1_0003 ## [5] 50 CACTCTGTGACCAGTGTGTGCGA…## [6] 50 CACTCTGTGACCAGTGTGCGCGA…50 CACTCTGTGATCAGTGTGCGCGA…AGTGCCCTGCCGGCATCGACTAT Hpl_8_0001 ## [8] 50 CACTCTGTGATCAGTGTGCGCGA…AGTGCCCTGCCGGCATCGACTAC Hpl_9_0001 ## [9] 50 CACTCTGTGATCAGTGTGCGCGA…AGTGCCCTGCCGGCACCGACTAC Hpl_10_0001 ## ## $nr ## [1] 464 62 39 27 37 33 54 248 20 ## ## $nm ## [1] 01 1 1 2 5 8 9 10

注意，单倍型出现在总数的2%以下，现在已经被删除。

filepath <执行(“extdata”、“Toy.GapsAndNs。fna"， package="QSutils")读取<- readDNAStringSet(filepath)读取

长度为100的DNAStringSet对象:## width seq names ## [1] 50 TGACGCGCACAGAGTGCTGCTA…Tgggttaccccgtcgtggncgc ## [2] 50 tgacgcgcacagagtgctgcta…Tgggttaccccgtcgtggtcgc ## [3] 50 tgacgcgcacagagtgctgcta…Tgggntaccccgtcgtggtcgc ## [4] 50 tgacgcgcacagagtgctgct -…Tgggttaccccgtcgtggtcgc ## [5] 50 tgacgcgcacagagtgctgcta…Tgggttaccccgtcgtggtcgc ## ... ... ...## [96] 50 tgacgcncacagagtgctgcta…tgggttaccccgtcgggtcgc ## [97] 50 tgacgcgcacagagtgctgcta…Tgggttaccccgtcgtggtcgc ## [98] 50 tgacgcgcacagagtgctgcta…Tgggttaccccgtcgtggtcgc ## [99] 50 tgacgcgcacagagtgctgcta…Tgggttaccccgtcgtggtcgc ## [100] 50 tgacgcgcacagagtgctgcta

(lstcollapse (reads) str <- DottedAlignment(lstcollapse $hseqs) data.frame(Hpl=str,nr= lstcollapse $nr)

# # Hpl nr TGACGCGCACAGAGTGCTGCTAAATGACTGGGTTACCCCGTCGTGGTCGC # # 61 # # 2  ......-...........................................2 # # 3  .....................-............................2 # # 4  ........................-.........................2 # # 5  ........................................-.........2 # # 6  ..........................................-.......2 # # 7……N ...........................................2 # # 8  -.......................................-.........1 # # 9  .-................................................1 # # 10  ..-...........................-...................1 # # 11  ....-.......-................-.................... 1 ## 12 ..........-....................................... 1 ## 13 ...............-.................................. 1 ## 14 ......................-.........N................. 1 ## 15 ........................N...........-.....N....... 1 ## 16 ..........................-....................... 1 ## 17 ...........................-..-................... 1 ## 18 ...............................-.................. 1 ## 19 ................................N................. 1 ## 20 ..............................................N... 1 ## 21 ................................................-. 1 ## 22 ................................................N. 1 ## 23 ..........................................N....... 1 ## 24 ..........................................N......N 1 ## 25 .........................................N........ 1 ## 26 ..................................N.N............. 1 ## 27 ...............................N.-................ 1 ## 28 ...........................N...................... 1 ## 29 .........................N........................ 1 ## 30 ......................N........................... 1 ## 31 ................N................................. 1 ## 32 .............N......................N............. 1 ## 33 ........N......................................... 1 ## 34 .....N............................................ 1

lstCorrected<-CorrectGapsAndNs(lstcollapse $hseqs[2:length(lstcollapse $hseqs)]， lstcollapse $hseqs[[1]]) #再次添加最丰富的单倍型。lstCorrected<- c(lstcollapse $hseqs[1]，lstCorrected) lstCorrected

长度为34的DNAStringSet对象:## width seq names ## [1] 50 TGACGCGCACAGAGTGCTGCTAA…Tgggttaccccgtcgtggtcgc 1 ## [2] 50 tgacgcgcacagagtgctgctaa…Tgggttaccccgtcgtggtcgc 2 ## [3] 50 tgacgcgcacagagtgctgctaa…Tgggttaccccgtcgtggtcgc 3 ## [4] 50 tgacgcgcacagagtgctgctaa…Tgggttaccccgtcgtggtcgc 4 ## [5] 50 tgacgcgcacagagtgctgctaa…Tgggttaccccgtcgtggtcgc 5 ## ... ... ...## [30] 50 tgacgcgcacagagtgctgctaa…Tgggttaccccgtcgtggtcgc 30 ## [31] 50 tgacgcgcacagagtgctgctaa…Tgggttaccccgtcgtggtcgc 31 ## [32] 50 tgacgcgcacagagtgctgctaa…Tgggttaccccgtcgtggtcgc 32 ## [33] 50 tgacgcgcacagagtgctgctaa…Tgggttaccccgtcgtggtcgc 33 ## [34] 50 tgacgcgcacagagtgctgctaa…TGGGTTACCCCGTCGTGGTCGC 34

lstRecollapsed < -Recollapse lstRecollapsed (lstCorrected, lstCollapsed nr)

## $nr ## [1] 100 ## ## $seqs ## DNAStringSet长度为1的对象:## width seq names ## [1] 50 TGACGCGCACAGAGTGCTGCTAA…CTGGGTTACCCCGTCGTGGTCGC 1

filepath <执行(“extdata”、“ToyData_FWReads。fna"， package="QSutils") lstFW <- ReadAmplSeqs(filepath,type="DNA") cat("Reads: "，sum(lstFW$nr)，"， Haplotypes: "，length(lstFW$nr)，"\n"，sep="")

##阅读:63736，单倍型:1153

filepath <执行(“extdata”、“ToyData_RVReads。fna"， package="QSutils") lstRV <- ReadAmplSeqs(filepath,type="DNA") cat("Reads: "，sum(lstRV$nr)，"， Haplotypes: "，length(lstRV$nr)，"\n"，sep="")

单倍型:1162

lstI < - IntersectStrandHpls (lstFW nr美元,lstFW工厂美元,美元lstRV nr, lstRV工厂美元)猫(“弗兰克-威廉姆斯和房车总写着:”,金额(美元lstFW nr) + (lstRV nr美元),“\ n”)猫(“弗兰克-威廉姆斯和房车读取上面的刺:“总和(lstI pFW美元)+ (lstI拟就美元),“\ n”)猫(“弗兰克-威廉姆斯单上面刺:“总和(lstFW nr /美元金额(lstFW nr美元)> 0.001),“\ n”)猫(“房车上面单刺:“总和(lstRV nr /美元金额(lstRV nr美元)> 0.001),“\ n”)猫(“\ n”)猫(“读入弗兰克-威廉姆斯:独特的单体型”,总和(lstI pFW美元[lstI拟就= = 0美元]),“\ n”)猫(“读入房车独特的单体型:”,总和(lstI拟就美元[lstI pFW = = 0美元]),“\ n”)猫猫(“\ n”)(“读取共同点:sum (lstI nr美元),“\ n”)猫(“单共同点:长度(lstI nr美元),“\ n”)

## FW和Rv的总读取:129256 ## FW和Rv的读取above thr: 99982 ## FW的单倍类型above thr: 35 ## Rv的单倍类型above thr: 35 ## ## FW唯一的单倍类型:874 ## Rv唯一的单倍类型:1047 ## ##常见的读取:98061 ##常见的单倍类型:1

4准物种数据探索

让我们加载一个玩具数据来演示不同的任务和函数。

filepath <执行(“extdata”、“ToyData_10_50_1000。fna"， package="QSutils") lst <- ReadAmplSeqs(filepath,type="DNA"

## $nr ## [1] 464 62 39 27 37 16 33 54 248 20 ## ## $hseqs ## DNAStringSet长度为10的对象:## width seq names ## [1] 50 caccctgaccagtgtgtgcga…gtgccccgttggcatgactat Hpl_0_0001|464|46.4 ## [2] 50 CACTCTGTGACCAGTGTGTGCGA…GTGCCCCGTTGGCATTGACTAT Hpl_1_0001|62|6.2 ## [3] 50 CACCCTGTGACCAGCGTGTGCGA…GTGCCCCGTTGGCATTGACTAT Hpl_1_0002|39|3.9 ## [4] 50 CACCCTGTGACCAGTGTGTGCGA…gtgccccgttggcatgactac Hpl_1_0003|27|2.7 ## [5] 50 CACTCTGTGACCAGTGTGTGCGA…GTGCCCCGTTAGCATTGACTAT Hpl_2_0001|37|3.7 ## [6] 50 CACTCGGTGACCAGTGTGCGTGA…GTGCCCCGTTGGCATTGACTAT Hpl_4_0001|16|1.6 ## [7] 50 CACTCTGTGACCAGTGTGCGCGA…GTGCCCTGTCGGCATTGACTAT Hpl_5_0001|33|3.3 ## [8] 50 CACTCTGTGATCAGTGTGCGCGA…GTGCCCTGCCGGCATCGACTAT Hpl_8_0001|54|5.4 ## [9] 50 CACTCTGTGATCAGTGTGCGCGA…GTGCCCTGCCGGCATCGACTAC Hpl_9_0001|248|24.8 ## [10] 50 CACTCTGTGATCAGTGTGCGCGA…GTGCCCTGCCGGCACCGACTAC Hpl_10_0001 20 | 2 |

的ConsSeq函数返回对齐产生的一致序列。该函数不考虑IUPAC歧义码，当出现平局时，将随机决定一致核苷酸。

ConsSeq (lst工厂美元)

[1] " cactctgaccagtgtgtgcgagtttggtgccccgttggcattgactat "

为了使构成准种的单倍型之间的差异形象化DottedAlignment函数返回一个字符串向量，每个字符串对应一个单倍型，其中一个点表示相对于主导单倍型的保守位点。

DottedAlignment (lst工厂美元)

# # Hpl_0_0001 | 464 | 46.4 # # # #“CACCCTGTGACCAGTGTGTGCGAGTTTGGTGCCCCGTTGGCATTGACTAT”Hpl_1_0001 | 62 | 6.2 # #“……T .............................................."39 # # Hpl_1_0002 | | 3.9  ## ".............. C ..................................."27 # # Hpl_1_0003 | | 2.7  ## "................................................. 37 C”# # Hpl_2_0001 | | 3.7 # #”……T .................................. 一个 ..........."16 # # Hpl_4_0001 | | 1.6 # #”……T.G ............ 让 ............................."33 # # Hpl_5_0001 | | 3.3 # #“……T .............. C……一个……T . . C ............"54 # # Hpl_8_0001 | | 5.4 # #“……T……………………T.CC…C……”# # Hpl_9_0001 | 248 | 24.8 # #“……T……………………T.CC…C…C”# # Hpl_10_0001 | 20 | 2 # # "…C T……………………T.CC CC…C”

的输出ReadAmplSeqs是否可以使用函数根据最丰富的单倍型的突变数进行排序SortByMutations，其中第一个单倍型是最相似的单倍型，最后一个单倍型是突变数量最多的单倍型。该函数还返回每个单倍型中突变数量的向量。

lstSorted < -SortByMutations (lst工厂美元,美元lst nr) lstSorted

长度为10的DNAStringSet对象:## width seq names ## [1] 50 CACCCTGTGACCAGTGTGTGCGA…gtgccccgttggcatgactat Hpl_0_0001 ## [2] 50 CACTCTGTGACCAGTGTGTGCGA…GTGCCCCGTTGGCATTGACTAT Hpl_1_0001 ## [3] 50 CACCCTGTGACCAGCGTGTGCGA…GTGCCCCGTTGGCATTGACTAT Hpl_1_0002 ## [4] 50 CACCCTGTGACCAGTGTGTGCGA…GTGCCCCGTTGGCATTGACTAC Hpl_1_0003 ## [5] 50 CACTCTGTGACCAGTGTGTGCGA…GTGCCCCGTTAGCATTGACTAT Hpl_2_0001 ##[6] 50仙人掌tgaccagtgtgcgtga…50 CACTCTGTGACCAGTGTGCGCGA…50 CACTCTGTGATCAGTGTGCGCGA…50 CACTCTGTGATCAGTGTGCGCGA…GTGCCCTGCCGGCATCGACTAC Hpl_9_0001 ## [10] 50 CACTCTGTGATCAGTGTGCGCGA…GTGCCCTGCCGGCACCGACTAC Hpl_10_0001 ## ## $nr ## [1] 464 62 39 27 37 16 33 54 248 20 ## ## $nm ## [1] 01 1 1 24 5 8 9 10

核苷酸或氨基酸在每个位置的频率可以用这个函数计算出来FreqMat．

FreqMat (lst工厂美元)

# # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17日18 19 20 21日22日23日24日25日26日# # 10 0 0 0 0 0 0 0 0 10 0 0 0 0 0 0 0 0 0 0 0 10 0 0 0 # # C 10 0 10 3 10 0 0 0 0 0 7 10 0 0 1 0 0 0 5 0 9 0 0 0 0 0 # # G 0 0 0 0 0 1 10 0 0 0 0 0 10 0 10 0 0 0 10 0 0 0 # # T 0 0 0 7 0 9 0 10 0 0 3 0 0 0 9 0 0 5 0 1 0 0 0 10 10 # # 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 # # 0 4 0 0 0 0 0 0 0 0 0 0 1 0 0 10 0 0 0 0 0 10 0 # # C 0 0 0 0 0 10 10 10 6 0 3 4 0 0 0 1 3 0 010 0 0 3 ## g 0 6 10 0 10 0 0 0 0 0 10 0 0 9 10 0 0 0 0 0 10 0 0 0 0 0 0 0 0 0 # t 10 0 0 10 0 0 0 0 0 0 4 0 7 6 0 0 0 0 0 9 7 0 0 0 0 10 0 7

为了在计算突变频率时考虑每个单倍型的丰度，单倍型丰度被传递给计算相同矩阵的函数，但具有丰度。

FreqMat (lst工厂美元,低水位体系域nr美元)

# # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 # # 0 1000 1000 0 0 0 0 0 0 0 0 0 0 0 0 # # 1000 C 1000 0 0 0 0 0 0 678 1000 1000 530 1000 0 0 39 0 # # G 0 0 0 0 0 1000 0 0 0 0 0 1000 1000 470 0 1000 # # T 0 0 0 0 322 984 1000 0 0 0 0 0 0 961 0 # # 17 18 19 20 21日22日23日24日25日26日27 28 29 30 31 32 # # 1000 0 0 0 0 0 0 0 0 0 0 355 0 0 0 0 # # 984 C 371 0 0 0 0 0 0 0 0 0 0 0 0 0 1000 # # G 0 1000 0 1000 1000 1000 0 0 0 0 0 1000 645 1000 0 0 # # 16 T 1000 0 629 0 0 0 0 0 0 1000 1000 1000 1000 0 0 # # 33 34 353.63.7 38 39 40 41 42 43 44 45 46 47 48 ## A 0 0 0 0 0 0 37 0 0 1000 0 0 0 1000 0 0 ## C 1000 1000 645 0 322 355 0 0 1000 0 20 322 0 0 1000 0 ## G 0 0 0 1000 0 0 963 1000 0 0 0 0 1000 0 0 0 ## T 0 0 355 0 678 645 0 0 0 0 980 678 0 0 0 1000 ## 49 50 ## A 1000 0 ## C 0 295 ## G 0 0 ## T 0 705

我们可能只对突变的体位感兴趣，所以MutsTbl得到的矩阵只报告每个位置突变的核苷酸或氨基酸的频率。

MutsTbl (lst工厂美元)

# # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17日18 19 20 21日22日23日24日25日26日27 28 29 # # 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 # # C 0 0 0 3 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # # G 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # # T 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 5 0 1 0 0 0 0 0 0 0 0 # # 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 # # 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 # # C 0 0 0 0 0 0 0 3 4 0 0 0 0 1 3 0 0 0 0 0 3 # # G 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # # T 0 0 0 0 0 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

正如可以做的那样FreqMat,MutsTbl单倍型的丰度可以传递给函数，得到一个具有丰度信息的突变表。

MutsTbl (lst工厂美元,低水位体系域nr美元)

# # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17日18 19 20 21日22日23日24日25日26日27 # # 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # # C 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 371 0 0 0 0 0 0 0 0 # # 16 G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # # T 0 0 0 322 470 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 16 0 0 0 0 0 0 # # 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 # # 355 0 0 0 0 0 0 0 0 0 0 37 0 0 0 0 0 0 0 0 0 0 0 # # C 0 0 0 0 0 0 0 0 0 322 355 0 0 0 0 322 0 0 0 00 295 ## g 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ## t 0 0 0 0 0 0 0 0 0 0 0 0 355 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

当序列特别大时，函数SummaryMuts计算一个表，显示排列中的多态性位置和观察到的每个核苷酸或氨基酸的频率。

SummaryMuts (lst工厂美元,lst nr美元= 0)

# # pos 530 C G T # # 4 4 0 0 470 # # 6 6 0 0 16 984 # # 11 11 0 678 0 322 # # 15 15 0 0 961 371年# # 19日19日0 0 0 984年629 # # 21 21日0 16 # # 28 28 355 645 0 0 # # 35 35 0 645 322 0 355 # # 37 37 0 0 0 678 # # 38 355 0 645 # # 39 39 37 963 0 # # 43 43 0 0 980 # # 44 44 0 322 0 678 # # 50 50 0 295 0 705

然后,PolyDist函数可以用更简单的方法得到多态位点的取代率。这个函数可以与丰度向量一起使用，也可以不使用。

PolyDist (lst工厂美元,低水位体系域nr美元)

## 4 6 11 15 19 21 28 35 37 38 39 43 44 ## 0.470 0.016 0.322 0.039 0.371 0.016 0.355 0.355 0.322 0.355 0.037 0.020 0.322 ## 50 ## 0.295

PolyDist (lst工厂美元)

## 4 6 11 15 19 21 28 35 37 38 39 43 44 50 ## 0.3 0.1 0.3 0.1 0.5 0.1 0.4 0.4 0.3 0.4 0.1 0.1 0.1 0.3 0.3 0.3

为了总结突变信息并计算每个突变的覆盖率，可以使用ReportVariants函数。该功能需要一个参考序列，在某些情况下可能是显性单倍型。

ReportVariants (lst工厂美元[2:长度(lst工厂美元)],lst工厂美元[[1]],lst nr美元)

# # 1 # # WT Pos Var浸C 4 T 879 # # 2 6 G 37 # # 3 C 11 C T 335 # # 4 62 T 19 C # # 388 # # 6 C 21 T 37 # # 7 G 28 351 # # 8 C 35 T 351 # # 9 T 37 C 335 T 38 C # # 351 # # 11 G 39 27 # # 12 T 43 C 248 # 335 # # 14 # 13 T 44 C T 50 C 341

另一种探索突变位置的方法是用每个位置的信息内容(IC)计算一个矩阵GetInfProfile．如果样品是DNA, IC最大值为2，而当样品是氨基酸时，IC最大值为4.32。

GetInfProfile (lst工厂美元,低水位体系域nr美元)

# # 1 2 3 4 5 6 7 8 # # 2.000000 2.000000 2.000000 1.002598 2.000000 1.881650 2.000000 2.000000 # # 9 10 11 12 13 14 15 16 # # 2.000000 2.000000 1.093457 - 2.000000 2.000000 - 2.000000 1.762312 - 2.000000 # # 17 18 19 20 21日22日23日24 # # 2.000000 2.000000 1.048563 - 2.000000 1.881650 - 2.000000 2.000000 - 2.000000 # # 25 26 27 28 29 30 31 32 # # 2.000000 2.000000 2.000000 - 1.061546 2.000000 - 2.000000 2.000000 - 2.000000 # # 33 34 35 36 37 38 39 40 # # 2.000000 2.000000 1.061546 2.000000 1.093457 1.061546 1.771636 2.000000## 41 42 43 44 45 46 47 48 ## 2.000000 2.000000 1.858559 1.093457 2.000000 2.000000 2.000000 2.000000 ## 49 50 ## 2.000000 1.124907

这可以被绘制出来:

dplot <- data.frame(IC=GetInfProfile(lst$hseqs,lst$nr)， pos=1:width(lst$hseqs)[1]) ggplot(dplot, aes(x=pos, y=IC)) + geom_point() + scale_x_continuous(minor_breaks =1: nrow(dplot)， breaks =1: nrow(dplot)) + theme(axis.text. frame)X = element_text(角度=45))

(#fig:plotic)对齐中每个位置的信息内容

5生物多样性指数对准物种复杂性的影响

随着准物种多样性指数的增加，这里有一个小插图:描述病毒准种也可在此包中使用。

6基因分型

另一个有趣的程序是对未知样本进行基因分型。为此，每个基因型都需要一套参考序列。建议至少有5个基因型特征良好的序列。这些集合应该是每个基因型的代表，并将提供基因型内方差的估计。函数DBrule，用于基因分型，考虑到目标单倍型和每个基因型之间的距离和基因型内的变异性。

第一步是装载目标单倍型用于基因分型ReadAmplSeqs，如下所示:

filepath <执行(“extdata”、“Unknown-Genotype。fna"， package="QSutils") lst2Geno <- ReadAmplSeqs(filepath,type="DNA") hseq <- lst2Geno$hseq[1] hseq . fna

长度为1的DNAStringSet对象:## width seq names ## [1] 285 CACGTCGCATGGAGACCACCGTG…100年18 GTTCAAGCCTCCAAGCTGTGCCT Hpl.0.0001 | |

然后使用相同的程序加载参考基因型序列。

filepath <执行(“extdata”、“GenotypeStandards_A-H。fas"， package="QSutils") lstRefs <- ReadAmplSeqs(filepath,type="DNA") RefSeqs <- lstRefs$hseq {cat("基因型引用序列的数量:\n") print(table(substr(names(RefSeqs)，1,1))}

基因型参考序列的数量:## ## A B C D E F G H ## 9 15 18 15 6 7 5 7

接下来，计算目标单倍型和参考单倍型之间的距离。参考单倍型之间的距离矩阵存储在dgrp中，在接下来的代码sniped中，而目标单倍型到参考序列的距离存储在向量d中DBrule函数根据两者计算出最可能的基因型，即从目标单倍型到参考文献的距离，以及同一基因型参考文献之间的距离。

dm <- as.matrix(DNA.dist(c(hseq,RefSeqs)，model="K80")) dgrp <- dm[-1，-1] d <- dm[1，-1] grp <- factor(substr(rownames(dgrp)，1,1)) hr <- as.integer(grp) dsc <- DBrule(dgrp,hr,d,levels(grp)) print(dsc)

## $Phi2 ## Phi2。Phi2。B Phi2. c Phi2. cD Phi2。E Phi2。F ## 0.0063243414 0.0057484677 0.0027527054 0.0006078818 0.0007949051 0.0188773955 ## Phi2。G Phi2.H DB.rule美元0.0448899802 - 0.0249913247 # # # # # # # #美元[1]4 # # # # # #型[1]“D”

目标序列被分类为基因型D，给出了最低的Phi平方值。

## R版本4.2.1(22-06-23)##平台:x86_64-pc-linux-gnu(64位)##运行在:Ubuntu 20.04.5 LTS ## ##矩阵产品:default ## BLAS: /home/biocbuild/bbs-3.16-bio /R/lib/libRblas. ##因此## LAPACK: /home/biocbuild/bbs-3.16-bio /R/lib/libRlapack。因此## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# [3] LC_TIME=en_GB LC_COLLATE= c# [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# [9] LC_ADDRESS=C LC_TELEPHONE= c# [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ## [8] base ## ##其他附加的包:## [1]QSutils_1.16.0 ggplott3 .3.6 ape_5.6-2 ## [4] Biostrings_2.66.0 GenomeInfoDb_1.34.0 XVector_0.38.0 ## [7] IRanges_2.32.0 S4Vectors_0.36.0 BiocGenerics_0.44.0 ## [10] BiocStyle_2.26.0 ## #### [1] tidyselect_1.2.0 xfun_0.34 bslib_0.4.0 ## [4] lattice_0.20-45 colorspace_2.0-3 vctrs_0.5.0 ## [7] generics_0.1.3 htmltools_0.5.3 yaml_2.3.6 ## [10] utf8_1.2.2 rlang_1.0.6 jquerylib_0.1.4 ## [13] pillar_1.8.1 withr_2.5.0 DBI_1.1.3 ## [13] glue_1.6.2 GenomeInfoDbData_1.2.9 lifecycle_1.0.3 ## [19] string_1 .4.1 zlibbioc_1.44.0 munsell_0.5.0 ## [25] gtable_0.3.1 psych_2.2.9 evaluate_0.17 ## [25] labeling_0.4.2 knitr_1. 1.0 ## [28] parallel_4.2.1 fansi_1.0.3 highr_0.9 ## [31]# [40] stringi_1.7.8 bookdown_0.29 dplyr_1.0.10 ## [43] grid_4.2.1 cli_3.4.1 tools_4.2.1 ## [46] bitops_1.0-7 magrittr_2.0.3 sass_0.4.2 ## [49] RCurl_1.98-1.9 tibble_3.1.8 crayon_1.5.2 ## [52] pkgconfig_2.0.3 assertthat_0.2.1 rmarkdown_2.17 ## [55] R6_2.5.1 nlme_1 .3 jsonlite_1.8.3 ## [37] farver_2.1.1 mnormt_2.1.1 digest_0.6.30 ## [40] stringi_1.7.8 bookdown_0.29 dplyr_1.0.10 ## [43] grid_4.2.1

#引用