1. 田间调查的一点小的思考
课题组研究的一大块内容是分子育种，正向遗传学方面的一些内容，主要的研究手段还是经典的 QTL + GWAS。这几天在处理大师兄的数据，大师兄主要用 GWAS 做产量性状的。
GWAS， 全基因组关联分析，主要是通过分析数量性状表型与基因型的关联统计分析挖掘候选基因。
处理数据的时候就发现了一些问题：
材料的系谱及名称。
我们在科学研究的过程中需要随时保持严谨，所用的材料来源，特性都应该有一个明确的记录，以便于后期重复使用以及鉴别自交子代的纯合程度。田间编号及材料管理
通常的我们在田间播种的时候会使用一套田间编号，便于播种管理，此编号通常有别于材料编号，因此需要留底，核对。我们在数据处理过程中发现，同一年同一地点同一个材料的两个重复的表型测量值差别较大，造成这种情况的原因推测有两种：a. 自交系材料杂合了; b. 田间编号有误。田间表型测量问题
表型测量是我们开展研究中比较重要的数据，一套好的表型数据可以运用到多方面的研究，因此特别需要注意表型值测量的准确。这里的准确主要指：a. 测量样本的均一；田间管理中发现材料的不同株系之间差别较为明显。b. 数据采集量足够，一般具有统计学意义的表型观测值通常 > 10 , 对于性状稳定的观测值，推荐不少于 5 个。田间表型调查额时机问题
由于调查的时机不同，如产量性状中的含水量，穗轴重，穗重等性状，随着时间的推移都会有一定程度的偏移，因此性状的调查一定要及时准确。一些细小的问题
性状调查前推荐制作好相应的性状调查表，数目一致，单位统一，测量方法一致，尽量减少人为引入的错误及误差。
2. 原始数据处理的一般办法：
原始数据调查完毕使用之前需要对原始数据进行初步的质控，主要检测由人工带入的一些错误数据，如录入时候小数点的错位等一些离群点。一般使用以下方法：
1. Tukey‘s test
E-value-min = Q1 - k(Q3-Q1)
E-value-max = Q3 + k(Q3-Q1)
k = 1.5 异常值
k = 3 极端异常值
对于 < 3 个观测值的测量样本，因其不具有统计学意义本方法及后续方法均不适用，因此建议田间调查观测值 一定要多，要准确。对于材料观测值较少的材料，进行人工干预(不推荐)，通过比较其在其他环境以及重复中的测量值确定其是否保留。缺失值以 NaN 表示。
perl 实现代码如下：
#!/usr/bin/perl -w
use strict;
use Data::Dumper;
use Statistics::Descriptive;
open IN, $ARGV[0] || die $!;
while(<IN>){
chomp;
next if /^$|^#/;
next if /sample/i;
my @a = split /\t/, $_;
my @tmp;
for (my $i=5;$i<@a;$i++){
if ($a[$i] !~ /NAN/i){
push @tmp, $a[$i];
}else{
next;
}
}
if (@tmp < 3 && @tmp > 0){
my $stat = Statistics::Descriptive::Full->new();
$stat->add_data(@tmp);
my $mean = $stat->mean();
print join("\t",@a[0,1,2,3,4], $mean)."\n";
}elsif(@tmp == 0){
print join("\t",@a[0,1,2,3,4], "NaN")."\n";
}else{
my $stat = Statistics::Descriptive::Full->new();
$stat->add_data(@tmp);
my $Q1 = $stat->quantile(1);
my $Q3 = $stat->quantile(3);
my $e_value_min = $Q1 - 3 * ($Q3 - $Q1);
my $e_value_max = $Q3 + 3 * ($Q3 - $Q1);
my @out;
foreach my $value (@tmp){
if ($value >= $e_value_min && $value <= $e_value_max ){
push @out, $value;
}else{
next;
}
}
my $out_stat = Statistics::Descriptive::Full->new();
$out_stat->add_data(@out);
my $out_mean = $out_stat->mean();
print join("\t",@a[0,1,2,3,4], $out_mean)."\n";
}
}
close;
__END__
Author : sxliulian2009@126,com
Date : 2018-01-11
这个程序可以剔除极端异常值，并最后返回观测值的均值，用于后续分析。
2. 格拉布斯准则法（Grubbs）
对于材料观测值数目 >10 的调查样本，我们可以采用该方法，保留 95% 置信区间内的值进行后续分析。
python 代码实现如下：
import numpy as np
from scipy.stats import t, zscore
def grubbs(X, test='two-tailed', alpha=0.05):
'''
Performs Grubbs' test for outliers recursively until the null hypothesis is
true.
Parameters
----------
X : ndarray
A numpy array to be tested for outliers.
test : str
Describes the types of outliers to look for. Can be 'min' (look for
small outliers), 'max' (look for large outliers), or 'two-tailed' (look
for both).
alpha : float
The significance level.
Returns
-------
X : ndarray
The original array with outliers removed.
outliers : ndarray
An array of outliers.
'''
Z = zscore(X, ddof=1) # Z-score
N = len(X) # number of samples
# calculate extreme index and the critical t value based on the test
if test == 'two-tailed':
extreme_ix = lambda Z: np.abs(Z).argmax()
t_crit = lambda N: t.isf(alpha / (2.*N), N-2)
elif test == 'max':
extreme_ix = lambda Z: Z.argmax()
t_crit = lambda N: t.isf(alpha / N, N-2)
elif test == 'min':
extreme_ix = lambda Z: Z.argmin()
t_crit = lambda N: t.isf(alpha / N, N-2)
else:
raise ValueError("Test must be 'min', 'max', or 'two-tailed'")
# compute the threshold
thresh = lambda N: (N - 1.) / np.sqrt(N) * \
np.sqrt(t_crit(N)**2 / (N - 2 + t_crit(N)**2))
# create array to store outliers
outliers = np.array([])
# loop throught the array and remove any outliers
while abs(Z[extreme_ix(Z)]) > thresh(N):
# update the outliers
outliers = np.r_[outliers, X[extreme_ix(Z)]]
# remove outlier from array
X = np.delete(X, extreme_ix(Z))
# repeat Z score
Z = zscore(X, ddof=1)
N = len(X)
return X, outliers
## PS: 代码转载，此代使用应人而异，使用前请适当修改。
异常值的检测方法还有很多，使用适合于自己数据的方法，才能在杂乱无章的数据海洋中打捞到自己的黄金。
总结
异常值的检测只是数据分析的开始，缺失值的插补更是统计分析中的一大神器。由于学识有限，掌握之日在再行补充。
田间调查工作强度高，要求严谨，是一件特别辛苦却大部分时候不被理解的工作，希望在田间调查的同学辛苦之余力求严谨，利人利己，这样科学研究的路才会越走越宽。
END
欢迎关注小刘哥
2018-01-17
