使用DMR進行PCA的主要目的是為了探究亞群中DMR是否存在大的差異，可以將不同亞群分開。
因為我們知道DMR在植物中存在三種類型CG、CHG和CHH，那么我們就需要對三種甲基化分別計算，以及總體的計算其PCA。

• 將馴化和改良分別得到的DMR進行合并
• 合并后從loci文件中計算每個個體在DMR的甲基化水平
• 合并甲基化水平并使用OSCA進行計算PCA

01 將馴化和改良分別得到的DMR進行合并

我們是兩次計算DMR，所以要計算總?cè)后wDMR的話，需要先對群體DMR進行合并，把DOM-DMR和IMP-DMR中overlap的地方避免重復計算，這個地方使用的是bedtools工具。

sort輸入文件

首先，以CG為例，我們把DOM-CG和IMP-CG進行重定向，然后要對重定向的新文件進行染色體和位置的排序，這樣處理會加快速度，需要的內(nèi)存也更少。并且如果文件沒有排序，bedtools merge就會報錯，所以sort是一定要的：

cat  DMR-CG.bed IMP-CG.bed > all_cg.bed

sort -k1,1 -k2,2n all_cg.bed > all_cg.sort.bed

合并之后可以使用bedtools 進行merge了，bedtools合并默認是有10個堿基的覆蓋，就在新的輸出文件中表示一次；也可以自己定義合并區(qū)間，比如你想對一定距離范圍內(nèi)的所有peak進行合并，那么需要使用-d 參數(shù)，例如合并500bp內(nèi)的peak，即bedtools merge -i all_cg.sort.bed -d 500

02 合并后從loci中計算每個個體的甲基化水平

因為我們是使用BatMeth2 進行甲基化的鑒定，這個軟件有一個比較人性的地方時默認給你輸出很多后需要用的甲基化格式文件，雖然也很占用儲存空間，但是對于一鍵式懶人來說還是比較不錯的，下面是其中一種輸出格式loci的結(jié)果：

loci.CG

每列分別是染色體，位置，方向，類型，甲基化的read，總C的read。
這里我是先將甲基化文件進行分不同染色體，并且計算甲基化的比率，笨拙的腳本如下：

import  os

#------------------------------------------------------------
#目的：獲取每個樣本DMR的甲基化程度
#
#01 將每個loci文件分成不同染色體的中間文件
#02 寫一個function 對每條染色體的甲基化程度進行統(tǒng)計
#03 合并所有每條染色體甲基化程度的中間文件，輸出最終結(jié)果，并刪除中間文件
#--------------------------------------------------------------

files = os.listdir()

for i in files:
    if  "loci" in i:
        #如果是 "loci"文件則進行打開
        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甲基化位點分為了12個不同染色體，接著我們分別對每個染色體進行計算，這樣會節(jié)省一些時間。

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])
    #設定好 染色體，起始 和 終止信息
    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)運行這個小腳本，這樣每條染色體DMR的甲基化都計算完畢。

合并甲基化水平并使用OSCA進行計算PCA

我們上兩步已經(jīng)得到了每條染色體的DMR甲基化水平，接著需要把數(shù)據(jù)合并，并整理成OSCA識別的格式。

for i in `cat sample.list`;
do
    cat ${i}_*ration.temp >> ${i}.CHG.meth &
done

接著把數(shù)據(jù)整理成這樣的格式即可：

OSCA輸入格式

接下來需要將多個meth文件進行合并，整理成最終osca可以識別的文件格式，具體腳本如下:

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軟件進行PCA的計算，計算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的格式后，使用其一個子命令進行PCA分析；而在PCA分析之前，需要計算樣品之間的親緣矩陣omics relationship matrix (ORM)，這個OSCA也給予了相關(guān)的命令，很簡單就可以完成：

osca --befile myprofile --make-orm --out myorm

得到的是三個同名的二進制文件，分別是bin,id和N.bin
最后終于到做PCA了，到這一步其實很簡單了，一行命令就可以解決：

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文件進行繪制PCA散點圖。

DMR-PCA