使用ADZE计算allelic richness (Ar) 和 private allele richness (pAr)

0.写在前面

软件官网:
https://web.stanford.edu/group/rosenberglab/adze.html
简单来说,等位基因丰度和私有等位基因丰度是描述品种遗传多样性的参数,大约类似于Ho与He。

sci论文.jpg

1.安装

简单的下载解压就行,提供linux和windows两种
tar -xzvf adze-1.0-XX.tar.gz

2.制作输入文件

需要两个输入文件,但本身并没有提供简单的创建输入文件的方法,因此我只能用简单暴力的方法尝试制作

1.datafile

目标格式.png

1.第一行是位点名称的列表,后面可以跟随一定数量的非数据行。然后每个个体的数据占用两行(对于单倍体或多倍体数据,每个个体的行数可能不同)。
2.在每一行的数据之前,有一定数量的列,用来表示个体的分类信息。图1中的文件在基因型数据之前有5列。在这个例子中,这些列分别代表(1)个体编号,(2)群体编号,(3)群体名称,(4)原产国家,和(5)原产地理区域。缺失数据的默认值是-9,但可以通过在参数文件中设置MISSING或者在命令行上使用-m选项来改变。
3.数据文件需要在参数文件中指定,并且必须与ADZE在同一个目录下;否则需要提供数据文件的完整路径。
4.等位基因的表示不限于数字:任何没有空白字符的字符串都可以用来表示一个等位基因。

1.质量控制后的vcf文件转fa格式

使用Northwest A&F University的Yudong Cai师兄提供的脚本
https://github.com/silvewheat/biotrans_cli/blob/master/vcf2fasta_hap.py

python3 vcf2fasta_hap.py -varfile test.vcf -outfile test.fa
fa文件.png

其中一组染色体一个序列,每个个体两条序列,因为不涉及单倍型问题,因此也不需要进行分型。对于不同的染色体,只是简单的首尾相连
注意:对于未知的点,标注的是“N”(第二行第一个)

2.拆分fa文件-每条序列一个文件

awk '/^>/ {s=substr($0,2)} {print > s".fa"}' test.fa

该脚本是将每条序列拆分,并以序列的表头文件命名,如:sample-1-1.fa,sample-1-2.fa,便于后续操作
建议将生成的fa文件移动至单独的文件夹,以免对总的fa文件造成影响

find . -name "*.fa" -exec sed -i '1d' {} \;

直接在原文件上,删除fa文件的第一行

3.多行文件转单行

find . -name "*.fa" -print0 | xargs -0 -I {} sh -c 'xargs < {} > {}.single'
单行序列.png

该命令是将换行符转换为空格,因此需要删除空格

find . -name "*.single" -exec sed -i 's/ //g' {} \;

尽管直接阅读显示多行,使用wc -l显示是单行。

4.每隔一个字符插入一个制表符

find . -name "*.single" -exec sed -i 's/./&\t/g' {} \;

5.将文件中N替换为"-9"

find . -name "*.fa.single" -exec sed -i 's/N/-9/g' {} \;

5.创建三个文件

1.序列编号文件
任意选择一个序列文件,按列编号

gawk 'NR==1 {n=NF} END {for (i=1;i<=n;i++) printf "rs%d\t", i; print ""}' sample.fa.single > num.txt

2.所有文件名的单列文件
按照sequence文件中的顺序合并fa.single文件

xargs cat < sequence > all.fa.single

3.最终datafile的前五列文件,与sequence一一对应

paste prefix all.fa.single >data
cat num.txt data > datafile

2.paramfile

文件名命名为:paramfile.txt
若要创建一个空模板ParamFile,请在默认ParamFile不存在且没有任何命令行参数的情况下运行Adze。 建议在模板的paramfile基础上设置参数

1.必须的参数

G 要计算的最大标准化样本量。 
DATA_LINES 数据文件中的数据行数(不包括包含位置名称的行) 
LOCI 位点数
NON DATA ROWS 在数据之前且应忽略的非数据行数。 此参数始终至少为1(位置名称的行)。 位置名称总是假定在第一行。 
NON_DATA_COLS 数据前面的非数据列数。 在图例中,非数据列数是5
GROUP_BY_COL 用于对个体进行分类的非数据列。 该值应小于NON_DATA_COLS且大于1
DATA_FILE 之前创建的输入文件的名字
R_OUT  等位基因丰富度输出文件的名称
P_OUT 私有等位基因丰富度输出文件的名称

2.组合参数

以下参数是计算组合群体的私有等位基因丰富度所必需的

COMB 默认设置为0,设置为1则表示需计算组合,必须设置以下几个参数
K RANGE 
原文这里有点难以理解,用newbing翻译了一下:如果数据分成了 A、B、C、D 四个组,那么从中选出两个组的组合有 AB、AC、AD、BC、BD、CD 六种。这段话要求所有的值(即组合的大小)必须大于等于 1,小于等于数据分组的个数(例如种群数)。如果组合的大小为 1,那么就相当于计算每个组的私有等位基因丰富度。可以用“-”字符来指定一个整数范围,例如 2-6 表示 2 到 6 之间的所有整数。整数和整数范围必须用逗号分隔,不能有重复的数字。
有效的规范示例有:2,2,5,7,2,4-6。这些示例分别表示要求 ADZE 计算所有两个组的组合的私有丰富度;所有两个组、五个组和七个组的组合的私有丰富度;以及所有两个组、四个组、五个组和六个组的组合的私有丰富度。
C_OUT 组合的私有等位基因丰富度输出文件的名称

3.高级参数

这些参数可以用来输出额外的信息,过滤掉那些有很多缺失数据的轨迹,或者指定缺失数据的替代方式。

MISSING 一个字符串,定义缺少的数据在数据文件中的表示方式。 默认设置是-9。 
TOLERANCE 这个参数是一个介于 0 和 1 之间的数值,用来指定在任何一个分组中,一个位点允许的最大缺失数据的比例。任何有至少一个分组的缺失数据比例超过这个数值的位点,都会被排除在计算之外。被排除的位点的名字会被写入到 deletedloci 文件中。这个参数的默认值是 1。
FULL R 这个参数是一个值为 0 或 1 的数值,默认为0。 为1则输出每个位点的等位基因丰度
FULL P 这个参数是一个值为 0 或 1 的数值,默认为0。 为1则输出每个位点的私有等位基因丰度
FULL_C 这个参数是一个值为 0 或 1 的数值,默认为0。 为1则输出组合的每个位点的私有等位基因丰度
PRINT PROGRESS 默认是1,表示在可能较长的计算期间输出进度报告。为0 则不输出

4.命令行参数

命令行标志为用户提供了从命令行输入信息的选项。 所有命令行参数都将覆盖ParamFile中指定的值。 如果指定了任何命令行参数,也必须给出ParamFile。

 -g MAX_G
 -d DATA_LINES
 -l LOCI
 -nr NON_DATA_ROWS
 -nc NON_DATA_COLS
 -s GROUP_BY_COL
 -f DATA_FILE
 -r R_OUT
 -p P_OUT
 -c COMB
 -k K_RANGE
 -o COUT
 -m MISSING
 -t TOLERANCE
 -tnocalc
 -fr FULL_R
 -fp FULL_P
 -fc FULL_C
 -pp PRINT_PROGRESS
使用方法:./adze-1.0 paramfile.txt -r richness_out
adze似乎需要完整路径

随后,运行./adze-1.0 paramfile.txt (adze需完整路径)
运行该命令所需内存约是文件大小的十倍

3.输出文件

1.R_OUT and P_OUT

P_OUT或R_OUT.png

这个文件包含了每个分组的等位基因丰富度(R OUT)或私有等位基因丰富度(P OUT)的汇总统计量。第一列是分组的名字(根据数据文件和 GROUP BYCOL 参数确定)。第二列是计算结果时使用的 g 值,最后三列是所有没有因为缺失数据而被过滤的位点的均值、方差和均值的标准误差。不同的分组之间用一个空行分隔。
对于P_OUT g值(也即样本量)是标准化的,所以不同组的g相同

2.C_OUT

C_OUT.png

这个文件包含了每个 k 个分组的组合的私有等位基因丰富度的汇总统计量。前 k 列是分组的名字(根据数据文件和 GROUP BYCOL 参数确定)。第 k+1 列是计算结果时使用的 g 值,最后三列是所有没有因为缺失数据而被过滤的位点的均值、方差和均值的标准误差。每个 k 个分组的组合之间用一个空行分隔。

3.FULL *

FULL_OUT.png

这个文件包含了等位基因丰富度(FULL_R)、私有等位基因丰富度(FULL_P)或每个 k 个分组的组合的私有等位基因丰富度(FULL_C)的汇总统计量和所有没有被过滤的位点的数据。第一行是所有用于计算的位点的名字。第一列或多列是分组的名字(对于 FULL_R 和 FULL_P)或多个分组的名字(对于 FULL_C)。第二列是计算结果时使用的 g 值。接下来的列是所有没有被过滤的位点的丰富度计算结果,最后三列是所有这些位点的均值、方差和均值的标准误差。每个分组或 k 个分组的组合之间用一个空行分隔。

4. _deletedloci

这个文件是在使用 TOLERANCE 参数来过滤有缺失数据的位点时生成的。这个文件的名字是 P OUT 文件的名字加上“ deletedloci”。它是一个报告,显示了由于指定的 TOLERANCE 而删除的位点的数量,以及删除的位点的名字的列表

5._summary

这个文件是在每次运行 ADZE 时生成的。这个文件的名字是 R OUT 文件的名字加上“ summary”。这是一个报告,显示了使用的参数和每个主要计算完成的时间(从执行开始算起)

6 例子

解压安装包后,其中以small为前缀的文件是示例

©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成,浏览时请结合常识与多方信息审慎甄别。
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

相关阅读更多精彩内容

友情链接更多精彩内容