Dsuite的使用 |
您所在的位置：网站首页 › f怎么占格子 › Dsuite的使用 |
|
Dsuite据说是相较于其他相似软件，占用的资源很少，这就给了我在自己服务器上操作的空间。 官网：https://github.com/millanek/Dsuite 官网上的教程：https://github.com/millanek/tutorials/tree/master/analysis_of_introgression_with_snp_data 1 安装要进行编译，必须拥有相当新的 GCC (>=4.9.0) 和 zlib 压缩库 (https://www.zlib.net），这也是为什么我装了新系统就去整了个gcc。 git clone https://github.com/millanek/Dsuite.git #直接克隆编译就行 cd Dsuite make #可选安装 cd utils python3 setup.py install --user --prefix= #=后面无任何符号，包括空格，用于画图 2 命令:Dtrios：计算D (ABBA-BABA)和f4-ratio统计的所有可能的种群/物种的三元 DtriosCombine：整合Dtrios在基因组区域运行的结果（如每条染色体） Dinvestigate：对D统计量显着提高的三元进行后续分析：计算f4统计信息，以及沿基因组窗口中的f_d和f_dM Fbranch（需进行上述可选安装）：计算与种群/物种相关的树的分支的D和f统计量 a) Dsuite Dtrios [OPTIONS] INPUT_FILE.vcf SETS.txt b) Dsuite Dquartets [OPTIONS] INPUT_FILE.vcf SETS.txt c) Dsuite Dinvestigate [OPTIONS] INPUT_FILE.vcf.gz SETS.txt test_trios.txt d) Dsuite Fbranch [OPTIONS] TREE_FILE.nwk FVALS_tree.txt 3 Dtrios 准备文件：vcf文件，接受bgzip和gzip压缩文件，可以包含Indel和多等位基因，但不会被计算（本文使用beagle文件是因为其他计算f4的软件都要求填充） sets.txt：两列文件，个体名称和种群名称 注意：必须要有一个Outgroup；对于不需要的种群，将种群名改为xxx即可屏蔽而不用提取vcf文件里的样本。 tree文件（可选，建议加上）：Newick 格式的树。这棵树应该有与物种/种群名称相对应的标签，与sets.txt文件中对应；分支长度可以存在但不会被使用（建议别加长度，可能出错）。需要注意的是，如果选择加上树文件，sets.txt文件里种群名必须是英文+数字的形式，可以有下划线。关于树文件的构建，教程放在末尾。 命令 ~/software/Dsuite/Dsuite/Build/Dsuite Dtrios test.beagle.vcf SETS.txt #可选命令： -k #默认=20，将数据集划分为Jackknife块的数量，整个数据集至少应该是20 -j #默认=NA，JJackknife块中snp的数量，如果使用了该命令，则替换-k -r，--region=start,length #只处理VCF文件的子集，如 --region=20001,10000 将处理 20001 到 30000 的变异位点 -t，--tree= TREE_FILE #一个带有newick格式的树的文件，该树指定种群/物种之间的关系。根据树排列的三组的D和f4-ratio值将输出在一个后缀为_tree.txt的文件中 -o，--out-prefix=OUT_FILE_PREFIX #指定输出文件的前缀，默认情况下，前缀取自SETS.txt文件的名称 -n，--run-name #run-name将包含在输出文件名的PREFIX后面 --no-f4-ratio #不计算f4-ratio -l NUMLINES #VCF输入的行数，如果通过unix管道读取VCF，则必须输入 -g，--use-genotype-probabilities # (optional) use probabilities (GP tag) or calculate them from likelihoods (GL or PL tags) using a Hardy-Weinberg prior，the probabilities are used to estimate allele frequencies in each population/species（原文贴在这，因为我看不懂） -p, --pool-seq=MIN_DEPTHVCF #包含pool-seq数据;也就是说，每个“个体”都是一个种群，然后从AD(等位基因深度)字段估计等位基因频率，例如MIN_DEPTH=5可能是合理的;当读取次数较少时，等位基因频率被设置为缺失 -c, --no-combine #不输出"_combine.txt"和"_combine_stderr.txt"文件，这些文件只有在DtriosCombine有用 --KS-test-for-homoplasy #检测强abba信息位点是否沿基因组聚集 并行程序此外，还提供了一个DtriosParallel命令，用法类似于Dtrios。允许同时输入多个vcf文件，该命令会调用DtriosCombine自动合并vcf文件，所有的vcf文件应包含相同的样本（如不同染色体的vcf文件） ./utils/DtriosParallel SETS.txt INPUT_FILE.vcf [INPUT_FILE.vcf ...] -j，-k，-t，-n，-l，-g，-p命令同上 -c，--no-combine #不运行DtriosCombine来获取单个综合结果文件 --cores CORES，#(default=CPU count) 设置进程数 --keep-intermediate #保留区域Dsuite Dtrios结果。 --logging_level {DEBUG,INFO,WARNING,ERROR,CRITICAL}, -v {DEBUG,INFO,WARNING,ERROR,CRITICAL} #最小级别的输出日志 --dsuite-path DSUITE_PATH #显式地将路径设置为Dsuite所在的目录。默认情况下，脚本将首先执行检查Dsuite是否可以从$PATH访问。如果不是，它将尝试在../Build/Dsuite中定位Dsuite。 --environment-setup ENVIRONMENT_SETUP, #应该运行该命令来为Dsuite设置环境。例如，'模块加载GCC'或'conda相较于Dtrios最大的优势应该是可以通过"--cores"设置进程数，其他大部分有用参数与上文类似。 结果文件head species_sets_no_geneflow_BBAA.txt P1 P2 P3 Dstatistic Z-score p-value f4-ratio BBAA ABBA BABA S01 S02 S00 0.00645161 0.228337 0.409692 7.6296e-05 40635 936 924 S01 S00 S03 0.0299321 1.08142 0.139754 0.000543936 28841 1724.75 1624.5 S01 S00 S04 0.0072971 0.308069 0.379015 0.00013279 28889.2 1691 1666.5 S01 S00 S05 0.00312175 0.133303 0.446977 5.6877e-05 28861.2 1687 1676.5 参考图.jpg 绘热图首先准备一个plot_order.txt文件，单列，包含每个种群名称。或者可以使用该命令： cut -f 2 sets.txt | uniq | head -n 20 > plot_order.txt其次准备两个ruby脚本文件 链接：https://github.com/millanek/tutorials/blob/master/analysis_of_introgression_with_snp_data/src/plot_d.rb https://github.com/millanek/tutorials/blob/master/analysis_of_introgression_with_snp_data/src/plot_f4ratio.rb ruby plot_d.rb SETS_BBAA.txt plot_order.txt 0.7 species_sets_no_geneflow_BBAA_D.svg ruby plot_f4ratio.rb SETS_BBAA.txt plot_order.txt 0.2 species_sets_no_geneflow_BBAA_f4ratio.svg 结果图.png该图是不考虑p1物种，列出p2和p3物种之间最大的D值和f4-radio值， 选择的种群有无基因流对结果有很大影响 4 DtriosCombine 整合上一步运行的结果（如不同染色体的运算结果） Dsuite DtriosCombine [OPTIONS] DminFile1 DminFile2 DminFile3 .... -o，--out-prefix=OUT_FILE_PREFIX #指定输出文件的前缀，默认情况下，前缀为“out” -t，--tree= TREE_FILE #一个带有newick格式的树的文件，该树指定种群/物种之间的关系。根据树排列的三组的D和f4-ratio值将输出在一个后缀为_tree.txt的文件中 -s , --subset=start,length #只处理VCF文件的子集，如 --subset=20001,10000 将处理 20001 到 30000 的变异位点 5 Dinvestigate 对D统计量显著升高的组进行后续分析 准备文件vcf文件和sets.txt：同上 test_trios：列出想要研究的组，可以多行，每行三列（意味着可以直接做想要的组而不必f4） #P1 #P2 #P3 mbuna deep Diplotaxodon对于这组，描述为：研究deep和Diplotaxodon组的混合信号，并和mbuna比较 命令 Dsuite Dinvestigate -w 50,25 test.beagle.vcf SETS.txt test_trios.txt 其中50 25是窗口和步长，对于动物来说，这个值似乎有点小了，建议自己摸索合适的参数 画图R脚本 plotInvestigateResults.R # read in the results with 50 SNP windows and a step of 25 SNPs bigStep |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |