关联分析LD计算的问题

说道关联分析，就要谈到连锁不平衡，连锁不平衡的衰减。 那么什么就是LD的衰减呢， 简单的说就是群体在经历了几百年或是更多年的历史重组之后在基因组上形成的重组块（block). 这些block 紧密的排列在基因组染色体的位置上。 有的地方block 比较大， 有的地方block比较小。 这个大小的意思就是延伸的长还是短。平均的延伸的长度决定了我们定位基因的精度。 延伸的越短，定位的精度越高。 延伸的越长，定位的精度越小。 我们画LD的衰减通常用tassel 软件， 或者是别的软件， 然后在进行做图计算。 大致分位一下几个步骤：

（1） 在TASSEL 中，计算LD时，会有三个选项，一个是full matrix， 一个是sliding window , 另一个是site by all . 这里我最熟悉的是full matrix 。 但是， 我们使用full matrix 的时候也有个缺点， tassel 不但计算了共线性的LD， 还计算了非共线性的LD， 大多的时候， 我们仅需要共线性的LD， 也就是一条染色体上两两位点之间的连锁不平衡大小。 当数据量比较大的时候， 最好按照染色体分开计算。 我喜欢的数据格式是hapmap.

(2)在进行计算之前，要进行位点的过滤。 选择 sites,MAF=0.05,remove indels,计算时选取为fullmatrix模式，这样计算的结果就是 同一条染色体上量位点的ld.

（3） 在R中首先过滤数据， 选取Dist_bp 和 R.2 那两列即可。 如果想用D` 值作为LD的度量的话，就可以选择D` 和Dist_bp 那两列。

（4） 根据Dist距离进行排序，有小到大。（此时，如果数据不是很大的话，就可以直接做图， 但通常当标记多的时候， 做的图，点很密，。有时看不到衰减的趋势）

（5） 按照窗口计算窗口内的R平方值的平均数。 这个窗口就是你打算按多少bp 来计算平均数， 就是每多少bp 计算一个平均数，沿着染色体依次计算。 可以是100, 1000,或者10000bp. 计算之后就可以按照窗口这个窗口设置的序列和平均的R 平方值来做图了。

以下是几个用到的R 代码， 计算窗口平均数的代码， 暂时还没写出来。 把分割hapmap 和进行清洗数据和整合的的代码附上。 计算窗口平均数的代码不是我写的。 所有无法附上。

# split hapmap file by each chromosomes
hap_split=function(hap){
hap_data<-read.csv(file=hap,sep="\t")
names(hap_data)
chr_names<-unique(grep("^\\d+", hap_data$chrom, value=T))
for ( i in 1: length(chr_names)){
  da_sub<-subset(hap_data, hap_data$chrom==chr_names[i])
  names(da_sub)[1]<-"rs#"
  names(da_sub)[6]<-"assembly#"
  filename=sprintf("chr%02d.hap.txt",i)
  write.table(da_sub,file=filename,  quote=FALSE, ,row.names=F,sep="\t")
}
}

## clearing the ld results from Tassel LD results , and all the results has put in one folder
cleardata<-function(file_adress){
  filename<-list.files(file_adress)
  for (i in 1: length(filename)){
    ld<-read.csv(file=filename[i], sep=",")
    ld_plot<-subset(ld,ld$pDiseq<0.01)  # select the significant pairwise ld under p<0.01
    ld_plot_data<-ld_plot[,c("Dist_bp","R.2")] # only select two cols containing Distance, RSquare
    names(ld_plot_data)<-c("Dist","R")
    ld_final<-ld_plot_data[order(ld_plot_data$Dist),]
    write_filename=sprintf("ld_chr%02d.hap.plot.txt",i)
    write.table(ld_final, file=write_filename, row.names=F,col.names = F, quote=F,sep="\t")
  }
}

combination<-function(file_adress){
  filename<-list.files(file_adress)
  com<-NULL
  for (i in 1:length(filename)){
  da<-read.table(file=filename[i],sep="\t")
  com<-rbind(com,da)
  }
  com<-com[order(com$V1),]
  write.table(com, file="ld_plot.txt", row.names=F,col.names = F, quote=F,sep="\t")
}