#####################################################################
# 基礎
#####################################################################

# 簡単な計算
1 + ( 2 - 3 ) * 4 / 5
sum( 1, 2, 3, 4, 5 )
help(sum)

# 変数と代入
x <- 5
x <- 2
x + 3

# データ型
a <- 12345
b <- "abcde"
c <- c( 1, 2, 3, 4, 5 )
d <- matrix( c( 1, 2, 3, 4, 5, NA ), ncol=2 )


#####################################################################
# RNA-Seq
#####################################################################

# パッケージインストール
#source("https://bioconductor.org/biocLite.R")
#biocLite("cummeRbund")


# データ読み込み ======================================================

## パッケージの呼び出し
library(cummeRbund) 

## サンプルデータ
file <- system.file("extdata", package="cummeRbund") #Rのファイル情報
cuff <- readCufflinks(dir=file)

## 自分のデータ
#cuff <- readCufflinks() # データ読み込み


# 2群間比較 ========================================================

## ステータス確認
replicates(cuff)

## 遺伝子レベルのデータの取り出し
cuffgene <- genes(cuff) 
## 散布図
csScatter(cuffgene, "iPS", "hESC")


# N群間比較 ========================================================

## クラスタリング
csDendro(cuffgene)

## 遺伝子IDベクトルの作成
data(sampleData) #サンプルデータセットの呼び出し
gene_vct <- sampleIDs #サンプルデータセットのピックアップ遺伝子
## 指定した遺伝子の情報の取り出し
cuffget <- getGenes(cuff,gene_vct)
## ヒートマップ
csHeatmap(cuffget)

## 遺伝子リストを自分で作成
# gene_vct <- c("XLOC_000004","XLOC_000005","XLOC_000008","XLOC_000009","XLOC_000011")


#####################################################################
# BS-Seq
#####################################################################

# パッケージインストール
#source("https://bioconductor.org/biocLite.R")
#biocLite("methylKit")

# データ読み込み ======================================================

## パッケージの呼び出し
library(methylKit) 

## 1sample
#file <- system.file("extdata","test.fastq_bismark.sorted.min.sam",package="methylKit") #Rのファイル情報
#file_list <- list(file)
#sample_list <- list("test") #ファイルに対応したサンプル名
#group_vct <- c(1) #ファイルに対応したグループナンバー
#bism <- processBismarkAln(file_list,sample_list,treatment= group_vct,assembly="hg18",read.context="CpG")

## 自分のデータ
#file_list <- list("test1.bam") #BAM（SAM）ファイル名のみ指定

## sample's
file1 <- system.file("extdata","test1.myCpG.txt",package="methylKit")
file2 <- system.file("extdata","test2.myCpG.txt",package="methylKit")
file3 <- system.file("extdata","control1.myCpG.txt", package="methylKit")
file4 <- system.file("extdata","control2.myCpG.txt", package="methylKit")
file_list <- list(file1,file2,file3,file4)
sample_list <- list("test1","test2","ctrl1","ctrl2") #ファイルに対応したサンプル名
group_vct <- c(1,1,0,0) #ファイルに対応したグループナンバー

bism <- methRead(file_list,sample_list,treatment= group_vct,assembly="hg18",context="CpG")

## 自分のデータ
#file_list <- list("test1.bam","test2.bam","control1.bam","control2.bam")


# メチル化比較 =======================================================

meth <- unite(bism) #統計用データの取り出し
diff <- calculateDiffMeth(meth)
methdiff <- getData(diff) #表データの取り出し
write.csv(methdiff,"sample.csv") #CSVファイルフォーマットで書き出し

# サンプル間の類似度検定 ==================================================

clusterSamples(meth) # サンプル間クラスタリングツリー
PCASamples(meth) # PCA解析