从TCGA下载的甲基化数据处理

上一篇已经写了从TCGA下载的甲基化数据的格式。
那么我想要知道某一个基因的甲基化情况要怎么办呢？
1、对下载数据放入一个文件夹
这一步在RNA数据的时候写过，用下面的脚本。

###用法为cmd perl 1.move_gzFileToFiles.pl
use strict;
use warnings;
use File::Copy;

my $newDir="data_in_one";
unless(-d $newDir)
{
    mkdir $newDir or die $!;
}

my @allFiles=glob("*");
foreach my $subDir(@allFiles)
{
    if((-d $subDir) && ($subDir ne $newDir))
    {
        opendir(SUB,"./$subDir") or die $!;
        while(my $file=readdir(SUB))
        {
            if($file=~/\.txt$/)
            {
                #`cp ./$subDir/$file ./$newDir`;
                copy("$subDir/$file","$newDir") or die "Copy failed: $!";
            }
        }
        close(SUB);
    }
}

2、筛选你所关注的基因的甲基化信息
这里自己写脚本

"""
作者：谢京合
时间：2020.10.28
功能：逐个搜索文件夹中的文件，并搜索匹配信息。
"""


import os


def main():
    """
    main function
    """
    path = 'D:/xiejinghe/data_in_one/'
    files = os.listdir(path)
    fresult = open('result.txt', 'w')
    # i_list = list[]

    for file in files :
        # print(file)
        fopen = open('./data_in_one/' + file, 'r')
        # print(file)
        flines = fopen.readlines()
        for i in flines:
            i = i.strip()
            i_list = i.split('\t')
            # if 'ITB1' in i:
            if 'ITB1;' in i and i_list[1] != 'NA':
                # i_list = i.split('\t')
                str_i = '\t'.join(i_list)
                # print('./test/' + file + '\t' + str_i)
                fresult.write(file + '\t' + str_i + '\n')
        fopen.close()
    fresult.close()

if __name__=='__main__':
    main()

3、将筛选得到的甲基化信息和对应基因的表达量进行合并和相关性分析
继续脚本处理，这里分为两个步骤。
第一步骤：因为一个样本对应很多甲基化位点信息，那么需要和基因表达量进行相关性分析之前，需要将每个样本的甲基化信息整合，用apply计算平均值或者中位值。

##第一步骤
setwd("E:/11.甲基化")

ITB1_methy <- read.csv("ITB1_methy.csv", header = T)
ITB1_methy_mean <- tapply(ITB1_methy$BetaValue,INDEX=ITB1_methy$id,FUN=mean)#这个可以求出来相同ID的均值。
ITB1_methy_median <- tapply(ITB1_methy$BetaValue,INDEX=ITB1_methy$id,FUN=median)#这个可以求出来相同ID的中位数。
write.csv(ITB1_methy_mean, file="ITB1_methy_mean.csv")
write.csv(ITB1_methy_median, file="ITB1_methy_median.csv")

第二步骤：将基因的表达信息和上一步处理完的甲基化信息合并，用merge进行合并。然后用cor进行相关性分析。

#第二步骤
library(plyr)
library(ggplot2)
library(export)
setwd("E:/11.甲基化")

ITB1_methy <- read.csv("ITB1_methy_mean.csv", header = T)
ITB1_expression <- read.csv("ITB1_expression.csv", header = T)

dataframe1 <- data.frame(ITB1_methy)
dataframe2 <- data.frame(ITB1_expression)
ITB1_frame1 <- merge(dataframe1,dataframe2,by="id",all=T)#表示全部要，但是两者之间差异的用NA填上。
ITB1_frame2 <- merge(dataframe1,dataframe2,by="id",all=F)#表示只需要两个数据集都有的部分，即交集。

ITB1_cor <- cor(ITB1_frame2[,2],ITB1_frame2[,3],method="pearson")#这个是-0.068
ITB1_cor <- cor(ITGB1_frame2[,2],ITB1_frame2[,3],method="spearman")#这个是-0.098


ggplot(ITB1_frame2,aes(x=betavale,y=log(value)))+ geom_point(size=1,shape=15)+geom_smooth(method=lm)
graph2ppt(file="ITGB1_cor.ppt", width=8, aspectr=1.5)

然后就得到了类似的结果。

image.png

完美的一天。