使用DMR進(jìn)行PCA的主要目的是為了探究亞群中DMR是否存在大的差異夭拌，可以將不同亞群分開口糕。
因?yàn)槲覀冎繢MR在植物中存在三種類型CG讨便、CHG和CHH，那么我們就需要對(duì)三種甲基化分別計(jì)算券盅，以及總體的計(jì)算其PCA。
- 將馴化和改良分別得到的DMR進(jìn)行合并
- 合并后從loci文件中計(jì)算每個(gè)個(gè)體在DMR的甲基化水平
- 合并甲基化水平并使用OSCA進(jìn)行計(jì)算PCA
01 將馴化和改良分別得到的DMR進(jìn)行合并
我們是兩次計(jì)算DMR膛檀，所以要計(jì)算總?cè)后wDMR的話锰镀，需要先對(duì)群體DMR進(jìn)行合并娘侍，把DOM-DMR和IMP-DMR中overlap的地方避免重復(fù)計(jì)算，這個(gè)地方使用的是bedtools工具泳炉。
sort輸入文件
首先憾筏，以CG為例，我們把DOM-CG和IMP-CG進(jìn)行重定向花鹅，然后要對(duì)重定向的新文件進(jìn)行染色體和位置的排序氧腰，這樣處理會(huì)加快速度，需要的內(nèi)存也更少翠胰。并且如果文件沒有排序容贝，bedtools merge就會(huì)報(bào)錯(cuò)，所以sort是一定要的：
cat DMR-CG.bed IMP-CG.bed > all_cg.bed
sort -k1,1 -k2,2n all_cg.bed > all_cg.sort.bed
合并之后可以使用bedtools 進(jìn)行merge了之景，bedtools合并默認(rèn)是有10個(gè)堿基的覆蓋斤富，就在新的輸出文件中表示一次；也可以自己定義合并區(qū)間，比如你想對(duì)一定距離范圍內(nèi)的所有peak進(jìn)行合并齐板，那么需要使用-d 參數(shù)厢蒜，例如合并500bp內(nèi)的peak，即bedtools merge -i all_cg.sort.bed -d 500
02 合并后從loci中計(jì)算每個(gè)個(gè)體的甲基化水平
因?yàn)槲覀兪鞘褂肂atMeth2 進(jìn)行甲基化的鑒定油额，這個(gè)軟件有一個(gè)比較人性的地方時(shí)默認(rèn)給你輸出很多后需要用的甲基化格式文件，雖然也很占用儲(chǔ)存空間刻帚，但是對(duì)于一鍵式懶人來說還是比較不錯(cuò)的潦嘶，下面是其中一種輸出格式loci的結(jié)果：
每列分別是染色體，位置崇众，方向掂僵，類型，甲基化的read顷歌，總C的read锰蓬。
這里我是先將甲基化文件進(jìn)行分不同染色體，并且計(jì)算甲基化的比率眯漩，笨拙的腳本如下：
import os
#------------------------------------------------------------
#目的：獲取每個(gè)樣本DMR的甲基化程度
#
#01 將每個(gè)loci文件分成不同染色體的中間文件
#02 寫一個(gè)function 對(duì)每條染色體的甲基化程度進(jìn)行統(tǒng)計(jì)
#03 合并所有每條染色體甲基化程度的中間文件芹扭，輸出最終結(jié)果，并刪除中間文件
#--------------------------------------------------------------
files = os.listdir()
for i in files:
if "loci" in i:
#如果是 "loci"文件則進(jìn)行打開
name = i.replace("res","")
for temp in range(0,13):
outname = name + str(temp)+"temp"
output = open(outname,"w")
for data in open(i).readlines():
my_data = data.strip().split("\t")
my_data[0] = my_data[0].replace("SL4.0ch0","")
my_data[0] = my_data[0].replace("SL4.0ch","")
if int(my_data[0]) == temp:
output.write(data)
這樣把總的CG甲基化位點(diǎn)分為了12個(gè)不同染色體赦抖，接著我們分別對(duì)每個(gè)染色體進(jìn)行計(jì)算舱卡，這樣會(huì)節(jié)省一些時(shí)間。
import sys
my_file = sys.argv[1]
outfile = my_file.replace("temp",".ration.temp")
out = open(outfile,"w")
for i in open("fin_cg_sort.bed").readlines():
my_list = i.strip().split("\t")
chr = int(my_list[0])
start = int(my_list[1])
end = int(my_list[2])
#設(shè)定好 染色體队萤，起始 和 終止信息
Me = 0
all_me = 0
for j in open(my_file).readlines():
data = j.replace("SL4.0ch0","")
data = data.replace("SL4.0ch","")
my_me = data.strip().split("\t")
if int(my_me[0]) == chr:
if int(my_me[1]) >= start and int(my_me[1]) <= end:
Me = Me + int(my_me[4])
all_me = all_me +int(my_me[5])
elif int(my_me[1]) > end:
if all_me == 0:
ratio = 0
out.write(str(chr)+'\t'+"cg"+'\t'+str(start)+'\t'+str(end)+'\t'+str(ratio)+'\n')
else:
ratio = Me / all_me
out.write(str(chr) + '\t' + "cg" + '\t' + str(start) + '\t' + str(end) + '\t' + str(ratio) + '\n')
break
else:
continue
循環(huán)運(yùn)行這個(gè)小腳本灼狰，這樣每條染色體DMR的甲基化都計(jì)算完畢。
合并甲基化水平并使用OSCA進(jìn)行計(jì)算PCA
我們上兩步已經(jīng)得到了每條染色體的DMR甲基化水平浮禾，接著需要把數(shù)據(jù)合并交胚，并整理成OSCA識(shí)別的格式份汗。
for i in `cat sample.list`;
do
cat ${i}_*ration.temp >> ${i}.CHG.meth &
done
接著把數(shù)據(jù)整理成這樣的格式即可：
接下來需要將多個(gè)meth文件進(jìn)行合并，整理成最終osca可以識(shí)別的文件格式蝴簇，具體腳本如下:
out=open("all_cg.site","w")
import os
all_file = open("all_cg.merge.bed").readlines()
all_line = len(open("all_cg.merge.bed").readlines())
files1 = os.listdir()
files = []
for i in files1:
if "meth" in i:
files.append(i)
out.write("Site")
for i in files:
if "meth" in i:
i= i.replace("_clean.CG.meth","")
out.write("\t"+i)
out.write("\n")
for i in range(0,all_line):
my_site = all_file[i].strip().split("\t")
site = "cg"+my_site[0]+"g"+my_site[1]
out.write(site)
for j in files:
if "meth" in j:
data1 =open(j).readlines()
for ii in data1:
this_site = ii.strip().split("\t")
if int(this_site[0])==int(my_site[0]) and int(this_site[2])==int(my_site[1]):
out.write("\t"+this_site[4])
out.write("\n")
接著使用OSCA軟件進(jìn)行PCA的計(jì)算杯活，計(jì)算PCA首先先把文件整理成BOD的格式：
osca_Linux --tefile myprofile.txt --methylation-m --make-bod --no-fid --out myprofile
#--efile reads a DNA methylation (or gene expression) data file in plain text format.
#--methylation-beta indicates methylation beta values in the file.
#--methylation-m indicates DNA methylation m values in the file.
#--make-bod saves DNA methylation (or gene expression) data in binary format.
#--out saves data ( or results) in a file.
#--no-fid indicates data without family ID.
整理成OSCA的格式后，使用其一個(gè)子命令進(jìn)行PCA分析熬词；而在PCA分析之前旁钧，需要計(jì)算樣品之間的親緣矩陣omics relationship matrix (ORM)，這個(gè)OSCA也給予了相關(guān)的命令互拾，很簡單就可以完成：
osca --befile myprofile --make-orm --out myorm
得到的是三個(gè)同名的二進(jìn)制文件歪今，分別是bin,id和N.bin
最后終于到做PCA了，到這一步其實(shí)很簡單了颜矿，一行命令就可以解決：
osca_Linux --orm myorm --pca 20 --out mypca
#--orm reads the ORM binary files.
#--pca conducts principal component analysis and saves the first n (default as 20) PCs.
最后用給出的eigenvec文件進(jìn)行繪制PCA散點(diǎn)圖寄猩。
