我是如何做生物统计作业的
上周生统课学习T检验，于是这周作业就是多重实验中如何运用T检验。刚好这题和上周介绍的数据管理有关系，所以放出来，求大家检查我的错误
题目
简单的说：对一个疾病的不同基因的两个样本进行T检验，然后修正p值。
原题如下：
Alzheimer’s Disease (AD) is a devastating neurodegenerative disorder affecting approximately 4 million people in the U.S. alone. AD is characterized by the presence of senile plaques and neurofibrillary tangles in cortical regions of the brain. These pathological markers are thought to be responsible for the massive cortical neurodegeneration and concomitant loss of memory, reasoning, and after aberrant behaviors that are seen in patients with AD.
Here we have a gene expression data from normal neurons and neurons containing neurofibrillary tangles of 14 mid-stage AD cases. The expression data can be found in , and description for samples in “GDS2795_sample.txt”. Use the information mentioned above to answer the following questions:
a) Use t-test to find significantly differential expression genes between normal and tangle neuron sample(p-value<0.05). Give the number of differential genes and the name of i th differential expressed gene ( i is the last two numbers of your student number, e.g. our student number is 201628010015316, then i=16, then the i th gene is “BC035139”).
Hint1: Notice that the two types of t-test with equal variance and unequal variance are different.
Hint2: “apply(data,1,function(x){…})” can apply function to every row in data more quickly than “for{}”. “names()” or “rownames()” can be used to extract names of differentially expressed genes
b) Adjust the p-values in question a) with both “” and “FDR” method to find differentially expressed genes in stringent way ( list the differentially expressed gene names and the adjusted p-value ).
Hint: you can do the adjustment according to the formula, or use “p.adjust()” instead.
注 题目所用的数据可以在我提供链接中下载：
解题思路
第一步都是先导入数据：
# 你需要把工作目录设置为数据存放点，或者指定绝对路径
expression_data <- read.table("GDS2795.txt")
sample_id <- read.table("/GDS2795_sample.txt")第二步简单审视数据，看看数据是怎么样子：
dim(expression_data)
#[1] 29849 14
head(expression_data)这样的工作对对另一个数据集也要进行，如果你用Unix类的系统，你可以直接在命令行下用我之前说的Linux数据分析里面提供的方法查看。
从简单的数据查看中，我们发现列名是病人，行名是各个基因。为了保持和sample_id的数据存放方式相同，我们需要进行一些数据清洗（也差不多是数据管理了）工作。
# 数据清洗
## 转置
t.expression_data <- as.data.frame(t(expression_data))
gene_number <- ncol(t.expression_data)
head(t.expression_data[1:5,1:5])
## 合并描述数据到表达数据中
t.expression_data$patient <- row.names(t.expression_data)
expr_data <- merge(t.expression_data,sample_id,by.x=c("patient"),by.y=c("V1"))
head(expr_data[1:5,1:5])
print(expr_data$patient)记住：每一步操作后都要检查一下数据的样子。bioinformatician的一大准则就是不要轻易相信你的数据结果，当然会提高二类错误(type II error)的概率了。
题目A
题目a要求是对每一个基因进行t.test，并且要根据方差是否相等，设置参数。显然不能简单的使用apply(x,FUN=t.test())。因此需要自己提供FUN
# 题目a：对每一个基因内两个样本进行假设性检验
# 考虑是否等方差
# 根据patient排序之后，得出normal和tangle依次出现
my.t.test <- function(x){
normal <- x[seq(from=1,to=14,by=2)]
tangle <- x[seq(from=2,to=14,by=2)]
p_value_f_test <-var.test(normal,tangle)$p.value
if (p_value_f_test > 0.05){
p_value_t_test <-t.test(normal,tangle,alternative = c("two.sided"),var.equal = TRUE)
}
else{
p_value_t_test <-t.test(normal,tangle,alternative = c("two.sided"),var.equal = FALSE)
}
return(p_value_f_test)
}
## 提取表达数据
expresion_gene_data <- as.matrix(expr_data[,3:ncol(expr_data)-1])
dim(expresion_gene_data)
expresion_gene_data[1:5,1:5]
results <- c(apply(expresion_gene_data,2,my.t.test))
hist(results)
rsults_order <- sort(results)
## 我的学号后两位是02,所以要选择第二名各个变量我使用我比较特有的错别字命名法，参考的时候前往注意，不要把我的错误拼写的代码也加进去。
题目B
其实做好第A题，就很简单了，第二题就是对结果进行修正而已。提示说了p.adjust()，我们就只需要简单的套函数就行了。
# 题目b:多重试验修正
# hint 使用p.adjust()
## 用法
#p.adjust(p, method = p.adjust.methods, n = length(p))
#p.adjust.methods
## c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
## "fdr", "none")
### bonferroni
p_adjusted_bonferroni <- p.adjust(results,method=c("bonferroni"))
table(p_adjusted_bonferroni < 0.05)
### fdr
p_adjusted_fdr <- p.adjust(results,method=c("fdr"))
table(p_adjusted_fdr < 0.05)结果我就得到一个比较显著性的基因，我也很是诧异，怀疑自己哪里错了。需要和TA的讨论一下
推荐阅读:
写在最后，希望看不到：
本来我希望使用dplyr和tidyr对数据处理，先把wide format 变成 long format，然后在进行聚合操作aggregate，但是无奈能力不行，还是继续 用了apply。
这就很尴尬了。