mecobalamin’s diary

人間万事塞翁が馬

Rで折れ線グラフ

Rで時系列データの折れ線グラフを作る
グラフを書くことよりも日付の処理が面倒だったのでそのメモ

作りたいのは気温と海水温の時系列グラフ
f:id:mecobalamin:20190320223529p:plain
気温と海水温の軸を左右に分けて、日付の軸は共通にする
折れ線グラフの赤が海水温で青が気温を示す

気温と海水温の時系列データは気象庁からダウンロードする
気温:気象庁、過去の気象データ
気象庁|過去の気象データ・ダウンロード
海水温:東京管区気象台
沿岸域の海面水温情報(関東・東海・北陸周辺) | 東京管区気象台

ファイルを見ると日付のフォーマットが違うことがわかる
今回に限らず、データによって日付の書き方が違うことはよくある
まず日付のフォーマットを揃える事から始める

気温のデータはcsv形式でダウンロードできて、
いろいろ情報が入っているが、
日付と気温以外のデータをexcelなどで削って使う
削ったデータをRで取り込んでみるとこんな感じ

data.air <- read.table("airtemp_2017-2019.csv", sep = ",")
>head(data.air)
1 ダウンロードした時刻:2019/03/20 12:34:56                                       
2                                                   羽田         羽田         羽田
3                                    年月日 平均気温(℃) 最高気温(℃) 最低気温(℃)
4                                  2017/1/1          8.6         12.5          3.7
5                                  2017/1/2          8.6         11.9          5.1
6                                  2017/1/3          8.9         13.6          3.7

これだと使いにくいのでヘッダーと何行スキップするかを指定する
また、ラベルを付け直して、一列目を日付のデータ形式に変換する
必要な列だけ取り出してデータの形式を変換する
気温のデータはファイルから読み込んだままだとfactor形式なので、
一度characterに変えてそれからnumericにする
手順についてググるといろいろ出てくる
factor型(2)_因子型(factor)から数値型(numeric)に戻すと数値が変わる? - 一所懸命に手抜きする

それと日付に関連してロケールを指定する必要があるらしい
なぜか?と聞かれると答えられないのだが。。。
コマンドはSys.setlocale()を使う

Sys.setlocale("LC_TIME", "C")
data.air <- read.table("airtemp_2017-2019.csv", header = TURE, skip = 2, sep = ",")
label.air <- c("date", "haneda_ave", "haneda_high", "haneda_low")
colnames(data.air) <- label.air
data.air$date <- as.POSIXct(as.Date(data.air$date, format = "%Y/%m/%d"))
data.air <- data.air[c("date", "haneda_ave")]
data.air[,1] <- as.Date(data.air[,1])
data.air[,2] <- as.numeric(as.character(data.air[,2]))

変換した結果がこれ

>head(data.air)
        date haneda_ave
1 2017-01-01        8.6
2 2017-01-02        8.6
3 2017-01-03        8.9
4 2017-01-04        9.5
5 2017-01-05        7.4
6 2017-01-06        5.8

海水温はテキストデータで手に入るが日付が
年・月・日で別の列になっていて、
ヘッダーと同じ文字列がデータの最後の行にも入っている

  yyyy mm dd areaNo. flag  Temp.
1 2017 01 01     704    R  25.16
2 2017 01 02     704    R  25.41
3 2017 01 03     704    R  25.58

ファイルを読み込んだら最後の行を取り除いて
年・月・日を一つにまとめたあとに日付のデータ形式に変換し、
必要な列だけ取り出す
こっちもfactorをcharacter -> numericの順に変換する

data.sea <- read.table("seatemp_2017-2019.txt", header = T, sep = ",")
data.sea <- data.sea[-nrow(data.sea),]
data.sea$date <- paste(data.sea[,1], data.sea[,2], data.sea[,3], sep = ".")
data.sea$date <- as.POSIXct(as.Date(data.sea$date, format = "%Y.%m.%d"))
data.sea <- data.sea[c("date", "Temp.")]
data.sea[,1] <- as.Date(data.sea[,1])
data.sea[,2] <- as.numeric(as.character(data.sea[,2]))

結果こうなる

        date Temp.
1 2017-01-01 25.16
2 2017-01-02 25.41
3 2017-01-03 25.58
4 2017-01-04 25.66
5 2017-01-05 25.66
6 2017-01-06 25.59

これを次の関数に投げるとグラフが出力される

PlotLineChart <- function(data.air, data.sea){
  # 日付の最小値と最大値を取得する
  Date.All <- as.Date(unique(append(data.air$date, data.sea$date)))
  Date.Range <- c(min(Date.All), max(Date.All))

  # マージンの調整、下、左、上、右の順
  par(mar = c(5, 5, 1, 5))

  # 1つ目のデータのプロット、軸を書かない
  # lは折れ線グラフ
  plot(data.air[,1], data.air[,2],
    xlab = "",
    ylab = "",
    xlim = Date.Range, # x軸の最小と最大
    ylim = c(-15, 35), # y軸の最小と最大
    type = "l",
    col = "gray",
    axes = FALSE
  )

  # 1つ目のデータのプロット、点を描く
  points(data.air[,1], data.air[,2], pch = 16, cex = 1, lwd = 1, col = "blue")

  # 軸ラベルを書く、4=右、3行目
  mtext("Air Temperature (Deg C)", side = 4, line = 3, cex = 1)
  axis(4, lwd = 1)

  # 重ね書きをする
  par(new = TRUE)

  # 軸だけ表示, type = "n"で指定する
  par(xaxt = "n")
  plot(0, 0,
    xlab = "Date (Year)",
    ylab = "Sea Temperature (Deg C)",
    type = "n",
    xlim = Date.Range, # x軸の最小と最大
    ylim = c(-15, 35), # y軸の最小と最大
    cex.lab = 1, # ラベルの文字の拡大率
    cex.axis = 1, # 軸の文字の拡大率
  )

  # x軸の表示
  par(xaxt = "s")
  axis.Date(1, at=seq(min(Date.All), max(Date.All), by = "years"), format = "%Y", cex.axis = 1)

  # 2つ目のデータのプロット、折れ線グラフで描く
  lines(data.sea[,1], data.sea[,2],
    cex = 1,
    col = "gray"
  )
  # 2つ目のデータのプロット、点を描く
  points(data.sea[,1], data.sea[,2],
    pch = 16,
    cex = 1,
    col = "red"
  )
}

グラフをファイルに出力するにはこうする

png.name <- sprintf("temperature.png")
png(png.name, width = 600, height = 400)
PlotLineChart(data.air, data.sea)
dev.off()

カラーパレットを作って月毎にマーカーの色を変える事もできる
まずは関数を作る

AddColorCode <- function(d, n, p){
  # カラーパレットを作る。並びを逆にしている
  Func.Col.Pal <- colorRampPalette(rev(RColorBrewer::brewer.pal(8, p)))
  # 12色のカラーコードを取り出して、輝度値をnにする
  Col.Pal <- Func.Col.Pal(12)
  for(i in 1:length(Col.Pal)){
    Col.Pal[i] <- paste(Col.Pal[i], n, sep = "")
  }

  tmp.name <- colnames(d)

  # 月の列を追加
  tmp.month <- as.numeric(substring(d$date, 6, 7))
  d <- cbind(d, tmp.month)

  # カラーコードの作成
  tmp.col <- c()
  for(i in 1:nrow(d)){
    tmp.col <- append(tmp.col, Col.Pal[tmp.month[i]])
  }
  d <- cbind(d, tmp.col)

  tmp.add.name <- c("month", "color")
  colnames(d) <- append(tmp.name, tmp.add.name)

  return(d)
}

データの前処理のときにAddColorCode()を呼び出して
データにカラーコードを追加する
3番目の引数はRColorBrewer::display.brewer.all()で表示される
カラーパレットの名前を使う
なおRColorBrewerはパッケージなので使う前にインストールが必要

data.air <- AddColorCode(data.air, 20, "Paired")
data.air <- data.air[c("date", "haneda_ave", "color")]
data.air[,1] <- as.Date(data.air[,1])
data.air[,2] <- as.numeric(as.character(data.air[,2]))
data.air[,3] <- as.character(data.air[,3])

data.sea <- AddColorCode(data.sea, 80, "Spectral")
data.sea <- data.sea[c("date", "Temp.", "color")]
data.sea[,1] <- as.Date(data.sea[,1])
data.sea[,2] <- as.numeric(as.character(data.sea[,2]))
data.sea[,3] <- as.character(data.sea[,3])

PlotLineChart()で色を指定するcolの引数を2ヶ所修正する
ここと

  points(data.air[,1], data.air[,2], pch = 16, cex = 1, lwd = 1, col = data.air$color)

ここ

  points(data.sea[,1], data.sea[,2],
    pch = 16,
    cex = 1,
    col = data.sea$color
  )

修正後にPlotLineChart()でグラフを表示させるとこんな結果になる
f:id:mecobalamin:20190320225337p:plain
見やすくなったとは言えないが例題ということで。。。。

JBrowseにplug-inを追加する

JBrowseにpluginを追加する

インストール方法はここを参考にした。
How do I install a plugin
JBrowse FAQ - GMOD

試しにSashimiPlotのプラグインをインストールする
まずpluginのダウンロードして展開
ディレクトリの名前をリネームするして
JBrowseのディレクトリにコピーする

SashimiPlot 0.1.3
Releases · elsiklab/sashimiplot · GitHub

$ wget https://github.com/elsiklab/sashimiplot/archive/0.1.3.tar.gz
$ tar -zxvf 0.1.3.tar.gz
$ mv sashimiplot-0.1.3/ /var/www/html/jbrowse/plugins/SashimiPlot/

次に/var/www/html/jbrowse/jbrowse.confにプラグインの名前を追加する
よくわからないけどディレクトリの名前と同じにする必要があった

[GENERAL]
plugins += SashimiPlot

/var/www/html/jbrowse/setup.shを実行する

$ ./setup.sh

最後にtrackList.jsonを編集する
sashimiplot/README.md at master · elsiklab/sashimiplot · GitHub

SashimiPlotで表示したいbamファイルを記入する

  {
      "label": "hogehoge_junctions",
      "storeClass": "JBrowse/Store/SeqFeature/BAM",
      "type": "SashimiPlot/View/Track/Sashimi",
      "urlTemplate": "hogehoge.bam",
      "chunkSizeLimit": 50000000,
      "useXS": true
  }

ブラウザで表示した結果がこれ
上からBigWigのヒストグラム、SashimiPlot、ch38のgff

f:id:mecobalamin:20190320145923p:plain

刺し身には見えないけれども。
Why is it called sashimi_plot?
What is sashimi_plot? — MISO documentation

scriptを止める

Atomが起動しなくなったので
再インストール

ここからダウンロード
Atom

インストールしたらそのまま使えた
パッケージもインストールされてた
設定ファイルは残っているみたいだ

起動しなくなった理由だけど
Rのスクリプトをscript pluginで実行したら
スペルミスがあって中断できなくなった
それでAtom終了させたら
Java Script errorが出て起動しなくなった

再インストールで復活したから良かったけど
scriptで実行しているプロセスの終了はctrl-qだったらしい

Command and shortcut reference

Command macOS Linux/Windows Notes
Script: Run cmd-i shift-ctrl-b If text is selected a "Selection Based" is used instead of a "File Based" run
Script: Run by Line Number shift-cmd-j shift-ctrl-j If text is selected the line number will be the last
Script: Run Options shift-cmd-i shift-ctrl-alt-o Runs the selection or whole file with the given options
Script: Run with profile shift-cmd-k shift-ctrl-alt-b Runs the selection or whole file with the specified profile
Script: Close View esc or ctrl-w esc Closes the script view window
Script: Kill Process ctrl-c ctrl-q Kills the current script process

beautiful soupで更新チェック

iPad mini 5が出るような噂
もしかしたらmini 4が安く出るんではないかと
整備品の値段をappleのサイトでいちいち調べてた

名前と値段のリストをとってくるスクリプトを書いてみた
PythonのBeautiful Soup4をつかう

まずiPadの正規整備品のサイトから
"refurbished-category-grid-no-js"のタグだけとってきて
その中に含まれる"li"のタグでfor文を回して
タグに含まれる製品名と価格を表示する

これだけだけど久々に書こうとしたら忘れてた

from urllib.request import urlopen
from bs4 import BeautifulSoup
import datetime

today = datetime.date.today()
todaydetail = datetime.datetime.today()

# date time
print('----------------------------------')
print("Today's refurbished iPad")
print(todaydetail)
print('----------------------------------')

url = "https://www.apple.com/jp/shop/refurbished/ipad"
html = urlopen(url).read()

bsObj = BeautifulSoup(html, 'lxml')

refurbishedipad = bsObj.find("div", {"class":"refurbished-category-grid-no-js"})

for i in refurbishedipad.find_all("li"):
    ipad_name = i.find("a").get_text()
    ipad_price = i.find("div", {"class":"as-price-currentprice as-producttile-currentprice"}).get_text()
    ipad_price = ipad_price.replace(" ", "").replace("\n", "")
    print(ipad_name + " " + ipad_price)

実行するとこんな感じ

----------------------------------
Today's refurbished iPad
2019-03-13 15:55:05.999841
----------------------------------
9.7インチiPad Pro Wi-Fi 32GB - スペースグレイ [整備済製品] ¥50,800(税別)
9.7インチiPad Pro Wi-Fi 32GB - ゴールド [整備済製品] ¥50,800(税別)
9.7インチiPad Pro Wi-Fi 32GB - ローズゴールド [整備済製品] ¥50,800(税別)
10.5インチiPad Pro Wi-Fi 256GB - シルバー [整備済製品] ¥73,800(税別)
10.5インチiPad Pro Wi-Fi 256GB - ゴールド [整備済製品] ¥73,800(税別)
10.5インチiPad Pro Wi-Fi 256GB - ローズゴールド [整備済製品] ¥73,800(税別)
12.9インチiPad Pro Wi-Fi 64GB - ゴールド [整備済製品] ¥73,800(税別)
12.9インチiPad Pro Wi-Fi 256GB - ゴールド [整備済製品] ¥87,800(税別)

atomで実行すると文字化けする
power shellで大丈夫だから
atomの問題っぽい

Atom-runnerの日本語文字化け解決方法 | げーまーこんぶの戯言

process.env.PYTHONIOENCODING = "utf-8";

追加したけど改善してない。。。

14 March 2019追記
PCの再起動とかしてたら文字化けしなくなった
PythonだけでなくRも大丈夫
エラーコードは相変わらず文字化けする

apt-getでファイルがダウンロードできない

apt-getでパッケージをインストールしようとしたら
ファイルがダウンロードできずにインストール出来なかったので
解決方法のメモ

欲しかったのはlibGL.so.1
opencv_createsamplesをインストールして実行したら
libGL.so.1がないってエラーが出た

パッケージの探し方もよくわからず。。。
ググるとパッケージ検索のサイトが見つかった
Ubuntu – Ubuntu Packages Search

"パッケージの内容を検索"
でキーワードに
"libGL.so.1"
を入れて
"キーワードに似た名前のファイルが含まれるパッケージ"
を探す
Ubuntu – Package Contents Search Results -- libGL.so.1

libgl1で良さそうなのでapt-getでインストール

$ apt-get install libgl1

だけど404エラーがでて必要なファイルがダウンロードできず。。。

解決方法がよくわからなかったけど、
apt-getのupdateで良かったらしい

$ sudo apt-get update

update後にapt-getを実行

$ sudo apt-get install libgl1

今回はエラーも出ず、
無事インストールできた

ググると出てきたのはインストール元を指定するファイルが
古くてダウンロードできないときの解決方法が多いみたいだった
apt-get updateが404で動かない場合もあるらしい
apt-get が 404 Not Found エラーで失敗する - eTuts+ Server Tutorial

メモ: リードカウントデータの正規化について

featureCountsでリードカウントしたあとに
複数のデータを比較する場合、
カウントの正規化が必要らしい

まだ良くわかってないので
参考にしているサイトのリンクを張っておく


ペアエンドのサンプルをリードカウントするときのコマンドがあった
featureCounts
RNA seqのリードカウント featureCounts - macでインフォマティクス

正規化の方法にはRLE正規化とTMM正規化があって、
求める正規化のファクターが異なる
これらの正規化方法を利用してるのがDESeq/DESeq2やedgeRで、
Rのライブラリーとして配布されている

DESeq/DESeq2はRLE正規化を使う
RLE正規化
RLE 正規化 | DESeq/DESeq2 に実装されている RNA-Seq カウントデータの正規化法

edgeRはTMM正規化を使う
TMM正規化
TMM 正規化 | edgeR で発現変動遺伝子を検出する際に利用する RNA-Seq データの正規化法

DESeq2を使った正規化でサンプルの処理をしている
サルマップ2018 (3) DESeq2による標準化と可視化まで
https://ncrna.jp/blog/item/388-deseq2-ggplot2

こっちはedgeR
edgeR - macでインフォマティクス

edgeRとDESeq2はbioconductorを使ってインストールする
Bioconductor - Home
パッケージマネージャーのBiocManagerをインストールしてから
ライブラリをインストールする

RはwindowsのRGuiを使っていて
versionは3.5.2

BiocManagerのインストール
Bioconductor - Install

edgeRのインストール
Bioconductor - edgeR

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("edgeR", version = "3.8")

DESeq2のインストール
Bioconductor - DESeq2

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2", version = "3.8")

とりあえずライブラリのインストールまで済んだ

hg38のfastaとgtf/gff3

GRCh38とhg38はヒトゲノムのシークエンスデータ
ほとんど一緒らしい。
前にも引用させてもらったけどここに説明がある。
GRCh38とhg38の違い(含むミトコンドリア)

GRCh38染色体名(ヘッダ行)が長く(gi|568815597|ref|NC_000001.11|)、hg38は簡潔(chr1)。

ヘッダ行に違いがあるのはやや問題で、
JBrowseで表示させたりfeatureCountsでリードカウントするときに
スムーズにできないことがあった

そこで簡潔な表示で揃えてfastaとgtf/gff3を作成する

まずはシークエンスデータのダウンロード
UCSC Genome Browser Downloads

リンクをたどるといくつかファイルが見つかる
chromosomeごとに分かれているデータを使う
Index of /goldenPath/hg38/bigZips

ファイルの説明はこちら

hg38.chromFa.tar.gz - The assembly sequence in one file per chromosome.
Repeats from RepeatMasker and Tandem Repeats Finder (with period
of 12 or less) are shown in lower case; non-repeating sequence is
shown in upper case.

ダウンロードと展開

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chromFa.tar.gz
$ gzip -d hg38.chromFa.tar.gz
$ tar -xvf hg38.chromFa.tar

展開するとパッチを含めて450ファイルぐらいある
パッチは取り扱いがよくわからないので
とりあえずchromosomeのファイルだけマージして使う

$ :> merge_chrom.fa
$ for i in {1..22}; do cat chr${i}.fa >> merge_chrom.fa; done
$ cat chrX.fa >> merge_chrom.fa
$ cat chrY.fa >> merge_chrom.fa

このfastaファイルをhisat2-build及びhisat2に使う
mecobalamin.hatenablog.com

次にhg38のgtfファイルを作成する
Gene annotation データを用意する(gtf形式) - Palmsonntagmorgen

NCBIからダウンロードできるgffファイルは詳しい表記のヘッダなので、
UCSCのサイトからgtfファイルをダウンロードしてgff3に変換する

Table Browser@UCSC
Table Browser

設定はこうした・・・ハズ。
clade: Mammal
genome: Human
assembry: Dec. 2013 (GRch38/hg38)
group: Gene and Gene Predictions
track: NCBI RefSeq
table: RefSeq All (ncbiRefSeq)
region: genome
output format: GTF - gene transfer format (limited)
output file: UCSC.hg38.gtf.gz
file type returned: gzip compressed

gtfファイルをgffreadをつかってgff3に変換する

$ gffread -E UCSC.hg38.gtf -o- > UCSC.hg38.gff3

このファイルはパッチで修正?された配列を含む
つまりパッチとchromosomeの間に重複がある
正しいやり方かよくわからないがパッチ側の重複した配列を除く
取り除かないと通らないコマンドがあった
bam2wigだったかな

$ cat UCSC.hg38.gff3 | awk '(a[$9]++ < 1000) && ($1 !~/[_alt]/) {print}' > rename_UCSC.hg38.gff3

exonやCDSで同じ名前のエントリーがあるので
1000個以上の重複で"_alt"を含む場合に除く
コードでは重複が1000個以下でかつ"_alt"を含まない行をファイルに出力している
1000個以上重複している配列はなさそうだったので、
使った感じでは取りこぼしはなさそう

今度はgff3をgtfに変換する

$  gffread rename_UCSC.hg38.gff3 -T -o- > rename_UCSC.hg38.gtf

gtfはfeatureCountsに、gff3はJBrowseに使う

ファイルの中身はこんな感じ
rename_UCSC.hg38.gtf

chr1    hg38_ncbiRefSeq exon    11874   12227   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    12613   12721   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    13221   14409   .       +       .       transcript_id "NR_046018.2"; gene_id "NR_046018.2";
chr1    hg38_ncbiRefSeq exon    14362   14829   .       -       .       transcript_id "NR_024540.1"; gene_id "NR_024540.1";
chr1    hg38_ncbiRefSeq exon    14970   15038   .       -       .       transcript_id "NR_024540.1"; gene_id "NR_024540.1";

rename_UCSC.hg38.gff3

# gffread -E UCSC.hg38.gtf -o-
# gffread v0.10.6
##gff-version 3
chr1    hg38_ncbiRefSeq transcript      11874   14409   .       +       .       ID=NR_046018.2;geneID=NR_046018.2
chr1    hg38_ncbiRefSeq exon    11874   12227   .       +       .       Parent=NR_046018.2
chr1    hg38_ncbiRefSeq exon    12613   12721   .       +       .       Parent=NR_046018.2
chr1    hg38_ncbiRefSeq exon    13221   14409   .       +       .       Parent=NR_046018.2
chr1    hg38_ncbiRefSeq transcript      14362   29370   .       -       .       ID=NR_024540.1;geneID=NR_024540.1
chr1    hg38_ncbiRefSeq exon    14362   14829   .       -       .       Parent=NR_024540.1
chr1    hg38_ncbiRefSeq exon    14970   15038   .       -       .       Parent=NR_024540.1