我是如何做生统作业的-答案修正版
2014级和2015级的生统考试发出去后，后面发现了一些小错误，虽然我明确强调只是提供思路参考，但是坑了那些对R不太懂，对统计也不太懂的小伙伴那就不太好了，于是我就重新把2年的题目整合在一起。这次应该齐全了。
题目一，基本无争议
研究者希望分析一个地区的Acute Keshan disease。健康人群和病人的phosphorus content结果放在了“Keshan_disease.txt”里，第一行是病人的结果，第二行是健康人群的结果。我们需要知道这两者是否有显著性的差异。
a. 在试验之前，你需要描述一下这两组之间的差异。比如说，平均值，画个箱图
b. 假设两组数据都来自于正态总体，检验一下均值是否有差异。你还需要检验两者的方差是否相同？
c. 假设没有第a题的假设，使用非参数方法去检验一下。
明确一下关键字： 描述性统计分析，方差检验，t检验，非参数检验。
请拿出《R语言实战》第二版，跟我一起翻书一起做。
第1小题请翻到第7章：基本统计分析。
首先是用summary了解他们的均值和四分位数
Keshan_disease_data <- read.csv("Keshan_disease.csv")
summary(Keshan_disease_data)
# patient healthy
# Min. :1.600 Min. :1.280
# 1st Qu.:4.228 1st Qu.:2.045
# Median :4.730 Median :3.675
# Mean :4.652 Mean :3.355
# 3rd Qu.:5.385 3rd Qu.:4.228
# Max. :6.530 Max. :5.780我还需要知道一下他们的标准差
apply(Keshan_disease_data,2,sd)
# patient healthy
#1.058033 1.298421各种统计函数，数学函数请翻到R语言实战的高级数据管理。这里划重点，我觉得apply是必考的内容，具体说明这里给出链接，不懂好好去看，看了还是觉得不懂的，你只能请我喝咖啡了。
画个箱图超级简单,请翻到R语言实战的基本图形：
boxplot(Keshan_disease_data)第2小题，首先要求我们检验方差是否相同, 用到是var.test。划重点，一定要先写你的假设，不然都不叫假设检验了。
H0: 健康的人样本的方差和病人样本的方差相同
Ha: 两组的方差不同
var.test(Keshan_disease_data$patient, Keshan_disease_data$healthy)
# F test to compare two variances
#data: Keshan_disease_data$patient and Keshan_disease_data$healthy
#F = 0.664, num df = 39, denom df = 39, p-value = 0.2055
#alternative hypothesis: true ratio of variances is not equal to 1
#95 percent confidence interval:
# 0.3511884 1.2554349
#sample estimates:
#ratio of variances
# 0.6639987由于P值> 0.05,我们有理由接受原假设(听了统计的文献课，我都不知道如何看待P值了)，因此方差相同。
然后我们就可以做t.test。不知道t.test怎么用怎么办，用help(test)，
Default S3 method:
t.test(x, y = NULL,
alternative = c(“two.sided”, “less”, “greater”),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95, …)
我们需要根据自己的要求修改参数，判断是否有差异，所以 alternative = c(“two.sided”), 方差相等，所以var.equal = TRUE。
还是记住，先把假设写了
H0: 两组的均值没有差异
Ha: 两组的均值存在差异
t.test(Keshan_disease_data$patient, Keshan_disease_data$healthy,
alternative = c("two.sided"), var.equal = TRUE)
# Two Sample t-test
# data: Keshan_disease_data$patient and Keshan_disease_data$healthy
# t = 4.8947, df = 78, p-value = 5.2e-06
# alternative hypothesis: true difference in means is not equal to 0
# 95 percent confidence interval:
# 0.7690202 1.8234798
# sample estimates:
# mean of x mean of y
# 4.65150 3.35525由于p=5.2e-06 小于0.01， 所以我们有把握拒绝原假设，认为两组的均值有差异（还是不知道如何看待P值）。
第3小题的非参数检验，建议用wilcoxon rank sum或 binomial test， 那我就用wilcoxon rank sum吧。请翻到基本统计分析的组间差异的非参数检验，方法就是wilcox.test，不知道怎么用，请用help(wilcox.test).
注意：写下你的假设，我这里就不写了。
wilcox.test(Keshan_disease_data$patient, Keshan_disease_data$healthy,
alternative = c("two.sided"))
# Wilcoxon rank sum test with continuity correction
# data: Keshan_disease_data$patient and Keshan_disease_data$healthy
# W = 1252.5, p-value = 1.358e-05
# alternative hypothesis: true location shift is not equal to 0
# Warning message:
# In wilcox.test.default(Keshan_disease_data$patient, Keshan_disease_data$healthy, :
# 无法精確計算带连结的p值发现了一个warning message， 你是不是有点慌，不要担心，只是警告而已，又不是error。如果不想看到warning添加一个exact = FALSE。
题目二： 基本无争议
“anemia.txt”给出了30个病人的anemia数据，我们希望用线性回归建立reticulyte和lymphocyte之间的模型。
注：lymphocyte是响应变量（也就是y了），请给出线性模型和p值。
关键字：线性模型
请翻到第8章：回归，用lm()拟合回归模型,对了记得看162页的函数哦。
anemia_data <- read.csv("anemia.csv")
anemia_fit <- lm(lymphocyte~reticulyte, data=anemia_data)
summary(anemia_fit)
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) -183.54 119.76 -1.533 0.137
# reticulyte 511.49 16.58 30.844 <2e-16 ***
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Residual standard error: 204.4 on 28 degrees of freedom
# Multiple R-squared: 0.9714, Adjusted R-squared: 0.9704
# F-statistic: 951.3 on 1 and 28 DF, p-value: < 2.2e-16所以p-value: < 2.2e-16,reticulyte的系数是511.49，那么模型就是
lymphocyte = 511.49 * reticulyte - -183.54
当然你还可以根据结果作图，看95%的置信区间，看皮尔森相关系数，看参数，对残差进行正态性检验等。
代码来自于计算所
plot(anemia_data$reticulyte, anemia_data$lymphocyte)
abline(anemia_fit)
#95% confidence interval for coefficients
confint(anemia_fit)
#pearson correlation coeefficient
cor(anemia_data$reticulyte, anemia_data$lymphocyte)
#residual
residual.v <- summary(anemia_fit)$residuals
#residual diagnosis with normality test
shapiro.test(summary(anemia_fit)$residuals)题目三： 争议在离群点
这次问题要用到datasets.txt， 数据来源于某个疾病的基因芯片。第一行是样品名（S1-S20）第一列是基因名（G1-G100)。数字都是表达量情况
a: G3基因在20个样本的密度图，（PDF哦），然后计算他的最小值，中位数，和方差。
b: 所有样本的基因分布的箱图。你需要检查一下离群点，然后删掉它重新来一遍。
这道题已经开始考验你的数据管理能力了，不太熟练的童鞋请到基本数据管理好好把代码敲一遍，当然我也写了专门的文章，这里是传送门。
A： 首先我们要用[]操作符提取数据
exp_data <- read.csv("datasets.csv", stringsAsFactors = F)
head(exp_data)
G3 <- exp_data[exp_data$X == "G3", ]
G3
# X S1 S2 S3 S4 S5 S6 S7 S8 S9
# 3 G3 9.881114 11.35708 10.27059 9.643856 10.37759 11.12393 10.55358 9.433481 9.793034
# S10 S11 S12 S13 S14 S15 S16 S17 S18
# 3 9.717762 10.13699 10.46827 9.600006 9.960798 10.02791 9.271463 10.46148 9.326767
# S19 S20
# 3 11.02746 11.05728结果G3是一个数据库， 我希望他是一个matrix或是一个vector，才可以让我方便作图，所以就要用到转置了。
G3.vector <- t(G3[2:length(G3)])那么问题来了，密度图，怎么画，急，在线求！！！请翻到120页核密度图
d <- density(G3.vector)
plot(d)min(G3.vector)
# [1] 9.271463
median(G3.vector)
# [1] 10.08245
var(G3.vector[,1])
# [1] 0.3880644PS: 核密度估计是用来估计随机变量概率密度函数的一种非参数方法。
B： 然后我们需要画箱图看看离群点（坑爹呢，输入法老提示利群店！！这广告也太明显吧！！）
boxplot(exp_data[,2:21])这个离群点还真不是一般的显著明显呀。先看这个数值大小
> exp_data$S7[exp_data$S7 > 80]
[1] 105
> exp_data$S8[exp_data$S8 > 80]
[1] 105估计是少打一个小数点把，没办法题目说删除，那我们也没有办法了，根据第72页做。
exp_data$S7[exp_data$S7 > 80] <- NA
exp_data$S8[exp_data$S8 > 80] <- NA
boxplot(exp_data[,2:21])争议点在于上图下面的一些点算不算离群点，还是就只有上面2点是离群点，提供一些其他方法：
# 直接在boxplot中排除离群点
boxplot（data，outline=FALSE）
# 先用boxplot.stats找到离群点，然后排除
boxplot.stats(ds3$S7)
which(ds3$S7 %in% boxplot.stats(ds3$S7)$out)
# 从boxplot中找到outlier，排除后作图
a<-boxplot(exp_data[,2:21])
a$out
na.omit(a$out) # 删除离群点
boxplot(a$stats)#取a里的数值重新作图题目四： 问题在于给定持家基因，要不要对proteinA标准化
蛋白A被人是肺癌病人预兆的预测信号。因此他们从三个医院收集了病人和健康的人数据，包括蛋白A的表达量，还有有beta-actin（持家基因），治疗后的存活时间等。数据在lung_cancer.txt。可以看instruction.pdf看详细细节。
a. 单因素ANOVA检验是否数据来源对蛋白的A表达量有影响
b. 一些研究者认为蛋白A表达水平是否异常与status没有关系。基于Alpha医院数据，我们可以否定这个结论么
提示： logistic回归或fisher.exact()
注意有一行代码是用来标准化的，lung_data$A.adujst <- lung_data$exp_A /lung_data$exp_actin， 新增一列标准化结果。
如果你不需要标准化，那么下面的A.adjust全部改为exp_A.
# 设置工作环境
setwd("c:/Users/Xu/Desktop/历年生统题目/2015生统考试/")
# 导入数据
lung_data <- read.table("lung_cancer.txt", header = TRUE)
# 浏览数据
head(lung_data)
# 标准化蛋白A，这样的标准化是不是有点简单粗暴，有没有更好的呢？
lung_data$A.adujst <- lung_data$exp_A /lung_data$exp_actin
# 之前在aov中因子化hospital，这里提前做
lung_data$hospital <- as.factor(lung_data$hospital)
# 看下结果
head(lung_data)
# 正态性检验
shapiro.test(lung_data$A.adujst[lung_data$hospital == 1])
shapiro.test(lung_data$A.adujst[lung_data$hospital == 2])
shapiro.test(lung_data$A.adujst[lung_data$hospital == 3])
# 方差齐性检验
bartlett.test(lung_data$A.adujst ~ lung_data$hospital)第A题： 请记得写假设哦
H0 : 不同医院来源的数据的蛋白A表达水平没有差异
Ha： 不同医院来源的数据的蛋白A表达水平存在差异
这题我建议你可以先去看下我写的方差分析（良心之作）。做方差分析脑子一定要画好一张表，如下
| 来源 | alpha | beta | gamma |
|---|---|---|---|
| x | 1 | 2 | 3 |
| … | … | … | … |
ANOVA分析，如果不要标准化, 那么把A.adujust 改为exp_A
exp_hospital.aov<-aov(A.adujst ~hospital, data = lung_data)
summary(exp_hospital.aov)
# Df Sum Sq Mean Sq F value Pr(>F)
# hospital 2 0.300 0.15023 7.284 0.000888 ***
# Residuals 197 4.063 0.02063
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1上次我到这里就结束了，上年的大神还做了多重比较, 发现第一个医院有问题。
TukeyHSD(exp_hospital.aov)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = A.adujst ~ hospital, data = lung_data)
$hospital
diff lwr upr p adj
2-1 -0.076735156 -0.13548034 -0.01798997 0.0065610
3-1 -0.078288694 -0.13703388 -0.01954351 0.0053863
3-2 -0.001553538 -0.06938663 0.06627955 0.9983886第2小题： 很简单，题目都说了fisher.test，翻到R语言实战143-144页。
fisher.test(lung_data$status[lung_data$hospital == 1], lung_data$abnormal[lung_data$hospital == 1] )
# Fisher's Exact Test for Count Data
# data: lung_data$status[lung_data$hospital == 1] and lung_data$abnormal[lung_data$hospital == 1]
# p-value = 0.008158
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
# 1.283453 9.026305
# sample estimates:
# odds ratio
# 3.338326或者是用glm来做，广义线性模型。
PA.alpha <- subset(PA, hospital == 1)
PA.alpha.fit <- glm(abnormal ~ status, data =PA.alpha, family = binomial())
summary(PA.alpha.fit)
# Coefficients:
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) -0.4595 0.3687 -1.246 0.21268
# status 1.2186 0.4502 2.707 0.00679 **题目五： 如何标准化数据，如何选择检验方法，如何计算power。
数据是2017考试版本的，也就是存在NA的那个版本。
标准化用的limma包，或者根据提示吧
检验方法：目前有wilcox.test或t.test
power: 后面放图
NASH是一种常见的肝病， 发现于肥胖和糖尿病有关。于是建立老鼠模型（可怜的小白鼠呀），研究adipose tissue和liver的关系。野生型雄性C57B1/6老鼠喂食LFD或者HFD，21周。然后根据liver histology和food fed分成4组。相信描述在GDS4013.txt, GDS4013_sample.txt和mehtylation.txt
A: 读取并标准化数据，看看有没有NA, 有的话，你如何处理
B： 找到HF>LF的上调基因，列出前20个，用啥检验呢？
C: 你有来自另一个研究的正常状态肝组织的每个基因的甲基化水平数据。你想知道差异表达基因的甲基化水平和其他正常状态的基因，请随机找DEG相同数量的数据，用t-test检验
D： 提供第三题的power.(温馨小提示：这是两个独立样本，当你计算绝对效应值的使用混池方差哦)
第一小题：
首先读取数据，然后看看有多少个缺失值。
nash <- read.table("c:/Users/Xu/Desktop/生统题目/考试时用到的数据/GDS4013.txt")
table(is.na(nash))
FALSE TRUE
481694 4发现有4个缺失值，不算多，可以直接把这几列的数据都删除（多重插补适合于许多数据缺失）。
which(is.na(nash), arr.ind = T)
row col
# PSPH 4 2
# TRDV2-2 26759 2
# PSENEN 10 3
# MC5R 26751 3
nash <- nash[c(-4,-26759,-10,-26751),]
table(is.na(nash))
# FALSE
# 481626新增： 我们也可以用na.omit()按行删除，就不用担心手动删容易出错了。
nash <- na.omit(nash)标准化应该会说明方法，这里我用limma包做。
library(limma)
nash.1 <- normalizeQuantiles(nash)第二题：我猜他的提示应该是说数据服从正态分布，要你用t.test作检验
首先计算var.test结果，放在最后一列
p.var <- apply(nash.1, 1, function(x) var.test(x[1:10], x[11:18])$p.value)
nash.2 <-cbind(nash.1, p.var)然后对每一行做t.test
# LF: x[1:10] HF: x[11:18]; HF>LF, LF <HF
# var.test result in x[19]
p.t <- apply(nash.2,1, function(x) t.test(x[1:10],x[11:18],
alternative=c("less"),
var.equal = x[19]>=0.05)$p.value)
table(p.t<0.05)前20个显著性差异的基因名。之前我放了20个显著性的p值，其实是手误了。
names(sort(p.t)[1:20])如果用非参数检验，估计大家都会做了。
第三题： 其实也没有多少问题
set.seed(1993)
DEG_name <- sample(names(p.t[p.t < 0.05]), 100)
non_DEG_name <- sample(names(p.t[p.t > 0.05]), 100)
methy_DEG <- methy[methy$V1 %in% DEG_name, ]
methy_non_DEG <- methy[methy$V1 %in% non_DEG_name, ]然后t.test, 我在里面附加了var.test.
t.test(methy_DEG$V2, methy_non_DEG$V2,
var.equal = (var.test(methy_DEG$V2, methy_non_DEG$V2)$p.value >= 0.05))
Two Sample t-test
data: methy_DEG$V2 and methy_non_DEG$V2
t = 4.7632, df = 198, p-value = 3.678e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.08116431 0.19585353
sample estimates:
mean of x mean of y
0.5340165 0.3955076第五题的power计算，估计也会继续考，这里更正一下我之前的代码
合并方差的算法是
n1 <- length(methy_DEG$V2)
var1 <- var(methy_DEG$V2)
n2 <- length(methy_non_DEG$V2)
var2 <- var(methy_non_DEG$V2)判断两个抽样的方差是否相同，选择不同的混池方法(如有疑问，马上和我说)
var.result <- var.test(methy_DEG$V2,methy_non_DEG$V2)$p.value
if (var.result > 0.05){
pool_var <- ((n1-1) * var1 + (n2-1) * var2) / (n1+n2-2)
}
else {
pool_var <- var1 / n1 + var2 /n2
}然后计算power。其中效应值 = (d1 -d2)/sd
library(pwr)
ds = (mean(methy_DEG$V2) - mean(methy_non_DEG$V2)) / sqrt(pool_var) # 根据混池计算效应值
powers <- pwr.t.test(n=100, d=ds, sig.level = .05,
type="two.sample", alternative = "two.sided")
powers题目六: 富集分析的计算方法
CD38 expression is an important prognostic marker in CLL with high levels of CD38 associated with shorter overall survival. In this study, we used gene expression profiling and protein analysis of highly purified cell-sorted CD38+ and CD38- chronic lymphocytic leukemia cells to elucidate a molecular basis for the association between CD38 expression and inferior clinical outcome.
Paired CD38+ and CD38- CLL cells derived from the same patient were used to perform the analysis.
(Please record your answer, R code, and R result)
Data: GDS2676.txt GDS2676_sample.txt
Please answer the following questions:
Read in data and do normalization. Draw a boxplot using all samples before and after normalization. (Hint: limma package, function ”normalizeQuantiles”)
Find differentially expressed genes (down-regulated, CD38- < CD38+) between CD38+ and CD38- disease samples, and provide top 20 down-regulated genes. (Hint: n is small, so please use non-parameter test)
Usually you are interested in the function indicated by differentially expressed genes, for which GO enrichment is a widely used method. In order to find whether the differentially expression genes (downregulated, p<=0.05) are enriched in “leukocyte activation during immune response” (GO term), please show a conclusive result using fisher exact test. Known genes annotated with this GO are listed in “GO_2_2.txt”
这是一道新题。题目说要我们用limma, 这是一个bioconductor包，不是基础R包。这就是给我们线索了，一旦给了非基础包，我们就可以知道会出哪一类的题目了。那么如何安装limma呢。
# install limma
source("https://bioconductor.org/biocLite.R")
biocLite("limma")我们先来学习一下limma
limmaUsersGuide(view=TRUE)然后你会看到146页的文档，你会惊呆了的。不过不用怕，我们就看quick start，其他部分都可以略过。大体上分为以下几步：
读取数据
预处理数据， 比如说标准化
拟合数据
查看结果
但是：根据题目的要求，我们连看都不用看了！！！！
基本上差异表达分析的数据读取之后都应该是如下这个样子
| \ | 样本A | 样本B | … |
|---|---|---|---|
| 基因A | x | x | … |
| … | … | … | … |
通过打开原文件，我发现体、提供的数据完全符合这个格式，根本不需要用到limma包专门的读取函数，直接用read.table就行了
第1小题：
data <- read.table("GDS2676.txt")
head(data)
boxplot(data)我们发现有很多离群点，而且这些离群点特别不一致，这就是为什么要做标准化的原因。题目说要用normalizeQuantiles, 我们可以看下这个函数的用法
library("limma")
?normalizeQuantiles
# Usage： normalizeQuantiles(A, ties=TRUE)
# A是一个matrix， 我们导入的数据符合
# ties：如果是TRUE，那么A的每一列都会仔细对待。好吧那就TRUE。然后标准化，画箱图：
library("limma")
?normalizeQuantiles
data.norm <- normalizeQuantiles(data)
boxplot(data.norm)大家的离群点特别的一致了，说明有些基因的确表达量很高。（这是的离群点就特别一一致，应该是不能排除把）
第2小题： 这时候要找差异表达的基因了，而且要用非参数的方法，而且还没有说用limma包，那就是我2014年题目中写的wilcoxon秩和检验的函数咯。
根据sample提供的数据，我们可以知道CDS-是1,3,5,7…. CDS+ 是2,4,6,8…，然后是CDS- < CDS+。还有他们来自于同一个细胞，所以应该是配对的。
这次换一种风格，用匿名函数和实体函数两种形式：
匿名函数
p.value<-apply(data2,1,function(x)wilcox.test(x[seq(1,11,2)],x[seq(2,12,2)],paired=TRUE,alternative = "less",exact=FALSE)$p.value)实体函数：
my.wilcox.test <- function(x){
x <- t(x)
CDS_sub <- as.numeric(x[seq(1,11,2)])
CDS_add <- as.numeric(x[seq(2,12,2)])
fit <- wilcox.test(CDS_sub, CDS_add, paired = TRUE， ternative = "less", exact = FALSE )
pv <- fit$p.value
return(pv)
}
pvalues.2 <- apply(data.norm, 1, my.wilcox.test)看看结果是否一致
all(pvalues == pvalues.2)应该返回TRUE， 返回FALSE的同学，请看下参数。
然后给出TOP 20下调基因，我认为这个TOP应该是要排序的，不然直接说给出前20个结果。
sort(pvalues)[1:20]第3小题： 这一题是要做富集分析，也就是填下表的过程。
| \ | 目标基因 | 全部基因-目标基因 |
|---|---|---|
| 注释 | A | C |
| 无注释 | B | D |
第一步算出目标基因的A和B。
table(DG.GENE %in% GO)
FALSE TRUE
607 10也就是说A=607，B=10。然后我们要从全部基因中把目标基因排除，然后计算C和D。
bg <- pvalues.2[! names(pvalues.2) %in% DG.GENE]
table(names(bg) %in% GO)
# FALSE TRUE
# 13295 170所以C是170，D是13295.
fisher <- matrix(c(10,607, 170, 13295), nrow=2)
fisher.test(fisher)
# Fisher's Exact Test for Count Data
# data: fisher
# p-value = 0.4593
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
# 0.6036833 2.4446923
# sample estimates:
# odds ratio
# 1.288372上课的PPT里面写了genome，还没有说背景是20300个基因，导致我理解偏差， 也怪我当时没好好听课。我觉得这次的结果应该是不会出错了。