前言
批次效应是测量结果中的一部分，它们因为实验条件的不同而具有不同的表现形式，并且与我们研究的变量没有关系。一般批次效应可能在下述情形中出现：
- 一个实验的不同部分在不同时间完成；
- 一个实验的不同部分由不同的人完成；
- 试剂用量不同、芯片不同、实验仪器不同；
- 将自己测的数据与从网上下载的数据混合使用；
- 不同建库策略，10X平台，Drop-seq, SMART2-seq
- 不同测序平台，BGI/Illumina
-
不同分析流程(甚至一个工具的多个版本，如salmon，CellRanger)
Leek JT, et al. Tackling the widespread and critical impact of batch effects in high-throughput data. Nat Rev Genet. 2010 Oct;11(10):733-9.
如果批次效应比较小还可以接受，如果批次效应很严重，就可能会和真实的生物学差异相混淆，让结果难以捉摸，我们需要对其进行校正。校正批次效应的目的就是：减少batch之间的差异，尽量让多个batch的数据相一致，这样下游分析就可以只考虑生物学差异因素。
已有的去批次方法：
- 芯片数据：RMA、MAS5、Z-score、分位数标准化等方法
- 高通量测序数据：R包sva的ComBat方法
- 单细胞转录组数据：Harmony、Seurat、LIGER等
Lazar C, et.al. Batch effect removal methods for microarray gene expression data integration: a survey. Brief Bioinform. 2013 Jul;14(4):469-90.
Tran HTN, et al. A benchmark of batch-effect correction methods for single-cell RNA sequencing data. Genome Biol. 2020 Jan 16;21(1):12
标准化（normalization）等方法只能缓解不能消除批次效应。
这些标化方法只考虑了少数影响批次效应的因素，通常并不能有效的去除批次效应，尤其是当处理跨实验室的样本数据时其标化效果更差，有些甚至可能加重技术性变异和扭曲真实的生物学信号
Wang D, et al. Extensive increase of microarray signals in cancers calls for novel normalization assumptions. Comput Biol Chem. 2011 Jun;35(3):126-30.
基因表达量的相对大小（秩序）关系
在正常组织样本中存在广泛且稳定的基因表达丰度排秩关系，可能反映了一种内在的生物学特性，为我们利用基因表达丰度排秩信息进行个体化差异表达分析提供了可靠的背景参照。
数据标准化对基因表达丰度排秩的影响
-
分别用Quantile 和Loess方法对数据进行了标化处理以及非标化处理，根据不同处理前后每个基因表达排秩的改变情况，绘制基因排秩变化的频数分布图来直观的反映标化对基因表达排秩的影响。
从标化原理上可以看出，Quantile方法（片间标化）不影响样本内基因表达排秩关系，Loess方法（片内标化）则会扰动样本内基因表达丰度的排秩关系。
1.基于基因表达量的秩序关系的差异基因分析方法RankComp
对于某一给定的疾病样本，我们利用每个基因相关的稳定基因对和相应的逆转基因对来判断它们是否发生差异表达。
对于基因i，如果在正常样本中它的排秩一致的高于基因j的排秩而在某一疾病样本中它们的排秩发生反转（即Gi > Gj → Gi < Gj），那么该基因对的逆转被认为是支持基因i在该疾病样本中发生差异下调。
相反的，如果在正常样本中它的排秩一致的低于基因j的排秩而在某一疾病样本中它们的排秩发生反转（即Gi < Gj → Gi > Gj），那么该基因对的逆转被认为是支持基因i在该疾病样本中发生差异上调。
为排除其它基因差异表达对基因i排秩的影响，假设上调和下调差异基因的比例是平衡的。
假定G表示基因i 相关的稳定基因对集合(i基因与其余基因），a和b分别表示两种不同排秩模式的稳定基因对（即Gi > Gj 和 Gi < Gj）；c和d分别表示两种不同反转模式的逆转基因对（即Gi > Gj → Gi < Gj和Gi < Gj → Gi > Gj）。
对于疾病样本k中某一特定的基因i，如果显著多的逆转基因对支持基因i上调，那么基因i 则被判断为差异上调；相反的,如果显著多的逆转基因对支持基因i 下调，那么基因i 则被判断为差异下调。
Fisher精确检验来判断支持该基因上调和下调的逆转基因对数目是否存在显著的区别，以此推断该基因是否发生显著的改变。在正常样本和疾病样本k中观察到Gi>Gj和Gi<Gj的频数分别是Obs11/Obs12=a/b和Obs21/Obs22=(a-c+d)/(b-d+c)。在零假设条件下（即支持该基因上调和下调的逆转基因对数目是相等的），我们期望观察到Obs11/Obs12=Obs21/Obs22。
基因i在样本k中差异方向的判定依赖两种不同反转模式的逆转基因对（即Gi >Gj→Gi<Gj 和Gi <Gj→Gi >Gj）：
如果c < d，差异上调；
如果c > d，差异下调；
否者，稳定表达。
R语言主要代码
kmp=matrix(c(reverse[,1],reverse[,2],reverse[,1]-reverse[,3]+reverse[,4], reverse[,2]-reverse[,4]+reverse[,3]),ncol=4)
GenePair_sig<-apply(kmp,1,function(x) fisher.test(matrix(x,ncol=2,byrow=T))$p.value)##fisher精确检验
2.基于基因表达量的秩序关系的差异基因分析RankcompV2
可以整合利用不同批次的小数据集进行分析
主要差别：对候选DEGs的过滤（过滤的目的是消除其他基因表达水平的改变对基因差异判断的影响）
RankComp中：对给定的基因G，其涉及的逆转基因对的伙伴基因是DEG，且差异方向与此逆转方向与此逆转基因对支持基因G的差异方向相反，则删除此基因对。
RankCompV2中：G配对的伙伴基因是候选DEG，则将此基因对从G特异性背景基因对中删除，重新统计正常和疾病样本中的基因对数目，做Fisher精确检验。迭代到两次临近迭代中的DEGs数目保持不变。
3.其他：基于秩序关系的通路分析方法individPath
3.其他：基于秩序关系的差异miRNA分析方法RankMiRNA
3.其他：基于秩序关系的差异lncRNA分析方法LncRIndiv
3.其他：基于秩序关系的血液细胞分子改变Ref-REO方法
3.其他：基于秩序关系的差异甲基化RMO方法
基因表达量秩序关系在非标准临床样本中的可重复性
病灶采样部位（不可控）
- 不同部位的癌细胞比例，差别大
-
活检甚至可能脱靶，癌旁
基因表达量秩序关系对癌组织取样部位（癌细胞比例）不敏感;
RNA降解、污染（难控）
- 样本预处理时，RNA极易发生部分降解；
-
FFPE样本核酸分子RIN约为2，远小于7。
基因表达量秩序关系对RNA降解不敏感 （FFPE 样本 VS FF样本）
基因表达量秩序关系在非标准临床样本中的可重复性
活检、CTC、单细胞微量问题(难控)
- 微量扩增偏倚
- 大量的缺失值
基因表达量秩序关系适用于微量样本, 如活检取样
基因表达量秩序关系的特点
- 平台（芯片、RNA-Seq）不敏感
- 无需数据片间标准化;
- 对癌组织取样部位（癌细胞比例）不敏感;
- RNA降解不敏感;
- 适用于微量样本;
Matlab实操
筛选稳定对
function stablepair=slect_stablepair(gid,exp,cutoff)
tic
len=length(gid)-1;
m=length(exp(1,:));
pair=nan(nchoosek(length(gid),2),3);
j=1;
for i = 1 : len
tic;
len-i
RE=exp(i*ones(len-i+1,1),:)-exp(i+1:end,:)>0;
ratio=sum(RE,2)./m;
pair(j:j+length(RE(:,1))-1,:)=[gid(i).*ones(len-i+1,1),gid(i+1:end,:),ratio];
j=j+length(RE(:,1));
toc;
end
stablepair=[pair(pair(:,3)>=cutoff,1:2);pair(pair(:,3)<=1-cutoff,[2,1])];
toc
end
# stable_pair_setA_90_lung=slec_stablepair(gid,'exp',0.9);
#%gid = GSE14520GPL3921geneid
#exp = GSE14520GPL3921
#cutoff = 0.9
#stable_pair_setB_90_lung=slec_stablepair('setB',0.9);
FDR校正
function [def_all,def_up,def_down]=fdr_adjust_bh_gqz(gene,p,index,fdr_cutoff)
m=length(p);
[p_sort,loc]=sort(p);
indexx(:,1)=(1:m).*(p_sort'<=(1:m)*fdr_cutoff/m);
n_sig=max(indexx);
if (n_sig)
sig_loc=loc(1:n_sig);
def_all=gene(sig_loc,:);
index=index(sig_loc);
def_up=def_all(index==1);
def_down=def_all(index==0);
else
def_all=[];
def_up=[];
def_down=[];
end
end
个体化差异分析
#[result,p_value,index]=individual_def_gene_sel5('con_pair_affy_illu_Agi_57148','gse27262_tumor',0.01);
function [result,p_value,index]=individual_def_gene_(stablepair,expE,fdr_cutoff)
#load(stable_pair)
#stable_pair=eval(stable_pair);
gidP=unique([stablepair(:,1);stablepair(:,2)]);
#load(exp_T)
gid=expE(:,1);
exp_tumor=expE(:,2:end);
for i = 1 : length(gidP)
length(gidP)-i
loc=find(stablepair(:,1)==gidP(i));
pair1=stablepair(loc,:);
a=size(pair1,1);
loc=find(stablepair(:,2)==gidP(i));
pair2=stablepair(loc,[2,1]);
b=size(pair2,1);
[~,loc1]=ismember(pair1,gid);
indexx1=exp_tumor(loc1(:,1),:)-exp_tumor(loc1(:,2),:)<0;
x=sum(indexx1,1);
[~,loc2]=ismember(pair2,gid);
indexx2=exp_tumor(loc2(:,1),:)-exp_tumor(loc2(:,2),:)>0;
y=sum(indexx2,1);
c=a-x+y;
d=b-y+x;
clear pair1 pair2 indexx1 indexx2 loc loc1 loc2
for j = 1 : size(exp_tumor,2)
table=[a,b;c(j),d(j)];
[~,p,~]=fishertest(table);
p_value(i,j)=p;
index(i,j)=(x(j))<y(j); #1 up; 0 down
end
end
for k =1 : size(exp_tumor,2)
time2=size(exp_tumor,2)-k
[def_all,def_up,def_down]=fdr_adjust_bh_gqz(gidP,p_value(:,k),index(:,k),fdr_cutoff);%0.05,1
result(k,1:3)={def_all,def_up,def_down};
end
end
