1前言
ROH参数的设置越看越复杂
《How to study runs of homozygosity using PLINK? A guide for analyzing medium density SNP data in livestock and pet species》这篇二区文章对各种参数的设置、有无过滤、Froh的计算做了详细的解释。但存在两个问题：20年发表的，只使用了芯片数据。
先说结论：如果只是想在文章中以三四句话的长度描述ROH和FROH,那么选取一篇同领域的前人文章，照猫画虎就足够了。
2.方法
主要使用的软件是最主流的PLINK
使用如下参数
plink \
--bfile ${input} \
--homozyg \
--homozyg-density 50 \ 一段ROH中每50kb必须有1个SNP
--homozyg-gap 100 \ 如果连续两个SNP的间隔大于100kb，那么就不能归为同一个ROH
--homozyg-kb 500 \ 只检测长度大于500kb的ROH
--homozyg-snp 50 \ 只检测长度超过50个SNP的ROH
--homozyg-window-het 1 \ ROH滑窗中可以允许有一个SNP位点为杂合
--homozyg-window-snp 50 \ 滑窗大小为50个SNP
--homozyg-window-threshold 0.05 \ 包含某个SNP的完全纯合滑窗的比例至少为5%
--homozyg-window-missing 5 \ 窗口缺失的snp数量
--homozyg-het 1 \ 允许ROH 片段里存在的最大杂合snp数量
--out ${output}
参数详解：
检测方法是扫描窗口法，分为三步
1.
扫描窗口首先由一个预设的snp数量定义，即上文的（ --homozyg-window-snp 50），最大的杂合SNP数量（–homozyg-window-het）和最大的缺失SNP数量（–homozyg-window-missing）则对这个窗口进行了限制。随后，根据这个定义的窗口扫描一个个体的基因组，并对每个snp评分，依据是每个snp出现在纯合窗口中的比例。（这里我的理解是步长为1的检测，对出现次数计数）
原文：The defined window stepwise scans an individual’s genome and scores for each SNP the proportion it appears in a homozygous window.
2.
使用snp评分的阈值识别纯合snp片段：
--homozyg-window-threshold 0.05：对于一个100个snp的窗口，一个snp至少要出现在5个纯合窗口，才会被识别为一个片段的一部分。
这种情况下，纯合窗口可能包含缺失或杂合的snp。
3.
最后，对这些纯合窗口设置额外的限制条件，首先检测一个片段中两个SNP之间的最大间隔（-homozyg-gap），以及最终ROH片段中允许的杂合SNP的最大数量（--homozyg-het）。不符合这两个条件的片段将被重新分割，这就会导致比滑动窗口更小的片段。最后，评估每个片段的最小SNP密度（以kb/SNP为单位）（--homozyg-density），以及最小长度和SNP数量（--homozyg-kb和--homozyg-snp）
注意：--homozyg-het参数在官网上是不存在的，我阅读的大部分文献也没有使用该参数。但是加入该参数，程序仍能正常运行）
3.Froh
对于常染色体的长度，一般是染色体上第一个snp到最后一个snp位置决定的，这在不同的实验中可能存在差异，此外，ROH片段长度还会受到上述参数的影响。
4.过滤参数
1.LD
在文中，作者建议通过改变参数（如snp数量和ROH长度）来避免假阳性ROH而不是进行LD过滤。
前人模拟研究建议应进行LD过滤，否则检测出来的共有片段很有可能是LD而不是ROH。但是本文发现随着LD修剪水平的提高，Froh迅速降低，LD过滤水平还会 影响到大片段的检测。
2.MAF
作者还建议不要在计算ROH前进行maf过滤，因为会大幅影响ROH片段长度。
MAF过滤在GWAS中被广泛使用，主要是因为GWAS非常强调单个基因型的精确度，但ROH中，单个snp的精度下降的影响很小。而且，实验结果表明，MAF过滤并不能改善ROH的检测。
5.检测参数
1.Minimal density(in kb/SNP)
结果表明，snp密度对ROH的影响很大
结果表明，PLINK默认的（50kb/snp）参数对大多数动物都是不合适的，可能会导致基因组覆盖的不完全。（这里作者使用的是芯片数据。对于重测序数据，则应该有很大不同，之前看到一篇文献使用的是10kb）。
对于不同的snp密度分布物种，作者建议使用（--homozyg-window-density）代替（--homozyg-density），强制算法检测每个窗口的snp密度而不是仅检查一个大片段。
2.Maximal gap
如图所示，即使是芯片数据，默认的1000kb也过大了，作者建议在保证最大化基因组覆盖的同时，最小化间隔。
3.Scanning window size and threshold
随着这两个参数的增加，ROH片段的数量会减少。
作者建议窗口的大小参数等于ROH的最小snp数量（L）。
对于阈值，作者建议使用如下公式计算
在另一篇文章中（https://doi.org/10.1186/s12711-020-00599-7），我找到了-window-threshold参数设置的计算方法
Nout<-2 #The number of outer SNPs you don't want to allow for ROH-analysis (calculation of windowthreshold)
l=log(0.05/(number_of_SNPs*number of animals))/log(1-mean_heterozygosity)#Lencz et al. 2007 Lencz, T., Lambert, C., DeRosse, P., Burdick, K. E., Morgan, T. V., Kane, J. M., ... & Malhotra, A. K. (2007). Runs of homozygosity reveal highly penetrant recessive loci in schizophrenia. Proceedings of the National Academy of Sciences, 104(50), 19942-19947.
这里的平均杂合度是使用plink的 "--hardy"计算，并对观测杂合度求平均值，忽略缺失值
l<-round(l)
if (is.na(l)){l<-50} #In a rare occasion l is set to NA, set l = 50 in this case
#Change window-snp to l
windowsnp<-l ##--homozyg-window-snp
#Calculate windowthreshold based on formula from Meyermans&Gorssen et al, 2020:
windowthreshold<-floor(1000*((Nout+1)/l))/1000
这里的脚本似乎是对不同群体分别计算的，我暂时不太能理解为什么不使用统一的参数
4.其他参数
其他作者使用的都是固定的参数
最小ROH长度为1mb，最大缺失snp为1，无杂合snp。
