###### tags: `ノート` `解析` `NCBI GEO` `R`

# 08_02： NCBI GEO解析　HeatMap作成 2

#### 2020/7/30 (月) , 在宅勤務

#### 【 作業場所＆保存場所 】
• //osc-fs_home/scratch/t-suzuki/disease_meth
• /home/s-maeda/osc-fs/08_NCBI_GEO_search
<br>

#### 再解析分も含めて、あらためてヒートマップ作成
きしまさん分のデータを `//osc-fs_home/scratch/t-suzuki/disease_meth`
に置いてもらったので、
きしまさん分のデータと、自分のデータのpeak.IQR.ratioをそれぞれで抽出し、結合してからヒートマップ作成した。

---

### 1：peak.IQR.ratioの抽出_1
まず自分の分から、
/home/s-maeda/osc-fs/08_NCBI_GEO_search に移動し確認。

・EPICディレクトリの数確認
`find ~/osc-fs/08_NCBI_GEO_search/GSE* -name EPIC | wc -l`
　98

・*mot_analysis_result.txtの数確認
`find ~/osc-fs/08_NCBI_GEO_search/GSE*/EPIC -name *mot_analysis_result.txt | wc -l`
　98

ディレクトリとtxtの数があっているので抽出

 #_抽出
for f in /home/s-maeda/osc-fs/08_NCBI_GEO_search/GSE*/EPIC/*mot_analysis_result.txt ; do
 awk '{print $14}' $f > $f.IQR
done

#_各列まとめる
paste /home/s-maeda/osc-fs/08_NCBI_GEO_search/GSE*/EPIC/*.IQR >> peak.IQR.ratio_pre.txt

#_max.peakヘッダー抽出
DIR="~/osc-fs/08_NCBI_GEO_search/GSE*/EPIC"
str=`echo ${DIR} | awk -F "/" '{ print $4 }'`
echo ${str} > peak.IQR.ratio_head.txt
sed -E 's/ +/\t/g' peak.IQR.ratio_head.txt > peak.IQR.ratio_head2.txt

#_ヘッダーと各peak.IQR.ratioあわせる
sed -i -e '1d' peak.IQR.ratio_pre.txt
cat peak.IQR.ratio_head.txt peak.IQR.ratio_pre.txt > peak.IQR.ratio_sm.txt

### 2：peak.IQR.ratioの抽出_2
次にきしまさんの分、
//osc-fs_home/scratch/t-suzuki/disease_meth に移動し確認。

・EPICディレクトリの数確認
`find //osc-fs_home/scratch/t-suzuki/disease_meth/GSE* -name EPIC | wc -l`
　34

・*mot_analysis_result.txtの数確認
`//osc-fs_home/scratch/t-suzuki/disease_meth/GSE*/EPIC -name *mot_analysis_result.txt | wc -l`
　528

*mot_analysis_result.txt の数が多い。
EPIC/ 以下のディレクトリにマージ前のサンプルがたくさんあるようだ。
こちらで消そうとしても、パーミッションの問題でできない。

とりあえず、自分のディレクトリにコピーして消した。

/home/s-maeda/osc-fs/08_NCBI_GEO_search/disease_meth へ移動し確認。

・EPICディレクトリの数確認
`find ~/osc-fs/08_NCBI_GEO_search/disease_meth/GSE* -name EPIC | wc -l`
　34

・*mot_analysis_result.txtの数確認
`find ~/osc-fs/08_NCBI_GEO_search/disease_meth/GSE*/EPIC -name *mot_analysis_result.txt | wc -l`
　34

たぶんOK。抽出。

 #_抽出
for f in /home/s-maeda/osc-fs/08_NCBI_GEO_search/disease_meth/GSE*/EPIC/*mot_analysis_result.txt ; do
 awk '{print $14}' $f > $f.IQR
done

#_各列まとめる
paste /home/s-maeda/osc-fs/08_NCBI_GEO_search/disease_meth/GSE*/EPIC/*.IQR >> peak.IQR.ratio_pre.txt

#_max.peakヘッダー抽出
DIR="~/osc-fs/08_NCBI_GEO_search/disease_meth/GSE*/EPIC"
str=`echo ${DIR} | awk -F "/" '{ print $5 }'`
echo ${str} > peak.IQR.ratio_head.txt
sed -E 's/ +/\t/g' peak.IQR.ratio_head.txt > peak.IQR.ratio_head2.txt

#_ヘッダーと各peak.IQR.ratioあわせる
sed -i -e '1d' peak.IQR.ratio_pre.txt
cat peak.IQR.ratio_head2.txt peak.IQR.ratio_pre.txt > peak.IQR.ratio_mk.txt

このディレクトリに、
`peak.IQR.ratio_sm.txt` と`motif_ID.txt` も移動した。確認もした。

---

### 3：ヒートマップの作成

#### ここからR

#データ読み込み
mk <- read.table("peak.IQR.ratio_mk.txt", header=TRUE, stringsAsFactors=TRUE)
sm <- read.table("peak.IQR.ratio_sm.txt", header=TRUE, stringsAsFactors=TRUE)
mo <- read.table("motif_ID.txt", sep = "\t", header = T)

#データ結合
mksm <- cbind(mk, sm)
data <- cbind(mo, mksm)

＃いちおう結合データをtxtにしておく
write.table (data, "peak.IQR.ratio_all.txt", sep ="\t", row.names=F, col.names=T, quote = F)


#パッケージ読み込み
library("ComplexHeatmap")    # ヒートマップを描く
library("circlize")　　　　# カラーパレットの色を細かく分ける
library("RColorBrewer")    # カラーパレット

#モチーフ名つける
rownames(data) <- unlist(strsplit(as.character(unlist(data["motif_ID"])), ".motif"))

#NAの行を削除
data_wona <- data[apply(data[,-1], 1, function(x){ !all(is.na(x))}),]

#色の設定
breaks1 <- seq(-2,2, length=500)

#breaks1に基づいて色を作製
mycols1 <- colorRamp2(breaks1,colorRampPalette(brewer.pal(9, "Blues"))(500))

#matrix変換
matrix1 <- as.matrix(data_wona)
matrix1 <- apply(matrix1,c(1,2), as.numeric)

#クラスタリング
df1 <- matrix1[,-1]
hc.rows <- hclust(dist(replace((df1), is.na(df1), 0)))
hc.cols <- hclust(dist(replace((t(df1)), is.na(t(df1)), 0)))

#ヒートマップ作成
ht1 <- Heatmap(df1,
                width = unit(6, "cm"),
                name = "peak/IQR",
                column_title = "Enriched motifs",
                column_title_gp = gpar(fontsize = 3),
                row_title = NA,
                row_names_gp = gpar(fontsize = 0.5),
                column_names_gp = gpar(fontsize = 0.5),
                col = mycols1,
                cluster_rows = hc.rows,
                cluster_columns = hc.cols,
                row_dend_width = unit(2, "cm"),
                na_col = "grey100",
                cell_fun = function(j, i, x, y, width, height, fill){
                    grid.rect(x = x, y = y, width = width, height = height, gp = gpar(col = "black", fill = NA, lty = 1, lwd = 0.2))
                }
 )

pdf("heatmap_20200730.pdf")
draw(ht1)

dev.off()