Palmsonntagmorgen

NGSデータを使った解析と、その周辺。

Linux ワンライナー覚書

Linux Genome NGS

スクリプトなどを用いずにLinuxのCUIコマンドのみで一行で行う処理（またはその処理を行うコマンド）のことをワンライナーと呼びます。ここではゲノム解析に使うワンライナーを順次追加していきます。

一般的なファイル操作

 $ ls workdir/            # workdirに含まれるファイルのチェック
 $ ls workdir/  | less  # ファイル数が多すぎる時はパイプしてlessで表示
 $ ls workdir/ | wc -l  # workdirに含まれるファイルの数を表示
 $ du -sh workdir/     # workdirディレクトリの容量チェック
 $ du -sh workdir/*   # workdirに含まれるファイルの容量チェック
 $ ls workdir/*.bed | wc -l  # workdirに含まれるBEDファイルの数チェック
 $ grep -x -f genelist1.txt genelist2.txt   # 2つのファイルの共通行を表示
 $ sed -i 's/aaa/bbb/g' sample.txt  # sample.txt 内の文字列aaaをbbbに置換
 $ sed -i 's/aaa/bbb/g' *.txt  # カレントディレクトリ内の全ての*.txt の文字列aaaをbbbに置換

BED, refFlatファイルなどのゲノムファイルの操作

BED

 # genome.faに含まれる染色体の表示
 $ grep \> genome.fa
 # BEDファイルをsorted BEDに変換
 $ sort -k1,1 -k2,2n peak.bed > peak.sort.bed
 # BED6ファイルをBED3に変換
 $ cut -f1-3 peak.bed6 > peak.bed3
 # peak.bedのピーク数をカウント（ヘッダが含まれる場合は-1）
 $ wc -l peak.bed 
 # peak.bedに含まれるピークを染色体別にカウント
 $ cut -f1 peak.bed | sort | uniq -c 
 # peak.bedからchr1のピークのみを取り出す
 $ cat peak.bed | awk 'if($1 == “chr1"){ print }' > peak.chr1.bed 
 # BEDファイル(タブ区切り)を２列目、１列目、３列目の順に表示
 $ cat peak.bed | awk -v 'OFS=\t' '{print $2, $1, $3}' 
 # カレントディレクトリ内にある全てのBEDファイルの染色体名 chrmt を chrM に置換
 $ sed -i 's/chrmt/chrM/gI' *.bed

refFlat

 # refFlatファイルをBED3形式に変換
 $ cat gene.refFlat | awk 'BEGIN{OFS = "\t"}{ print($3, $5, $6)}' > gene.bed 
 # refFlatファイルをBED6形式に変換
 $ cat gene.refFlat | awk 'BEGIN{OFS = "\t"}{ print($3, $5, $6, $1, ".", $4)}' > gene.bed
 # refFlatファイルからForward strandの遺伝子のみを取り出す
 $ cat refFlat.txt | awk 'if($4 ==+){ print }' > refFlat.forward.txt
 # refFlatファイルからReverse strandかつchr2の遺伝子のみを取り出す (複数条件指定)
 $ cat refFlat.txt | awk '{if($3 == “chr2” && $4 ==“-”){ print }}'  > refFlat.reverse.chr2.txt
 # refFlatファイルに含まれる遺伝子のうち最大のexon数を表示
 $ cat refFlat.txt | grep -v exon | awk 'BEGIN{ a=0 }{ if(a<$9) a=$9} END{ print a}'
 # refFlatファイルに含まれる遺伝子の平均exon数を表示
 $ cat refFlat.txt | grep -v exon | awk 'BEGIN{ a=0 }{ a+=$9 } END{ print a/NR }'
 # 平均exon数を整数で表示したい場合
 $ cat refFlat.txt | grep -v exon | awk 'BEGIN{ a=0 }{ a+=$9 } END{ printf "%d\n", a/NR }'

中級者以上向け

 # ファイル名にハイフンを含むBEDファイルのハイフンをアンダーバーに一括変換
 $ for file in `ls *-*.bed`; do mv $file ${file/-/_}; done
 # カレントディレクトリのsraファイルを全てfastq-dump
 $ find . -name '*.sra' -exec fastq-dump {} \; 
 # xargsを使って上と同じ処理
 $ ls *.sra | xargs fastq-dump
 # カレントディレクトリ以下にある全ての.txtファイルに対してaaaからbbbへの置換を実行（元のファイルは*.bakとして保存される）
 $ find ./ -name '*.txt' | xargs sed -i.bak "s/aaa/bbb/g"

解析ツールを利用したワンライナー

# ヒトゲノム (hg38)の染色体長を表すBEDファイルを作成
$ fetchChromSizes hg38 | awk -vOFS="\t" '($1!~/_/){ print $1, "0", $2 }' \
   > chrlen.hg38.bed
# merged siteが各ピークファイルといくつ重複しているかを調べる
$ intersectBed -wa -a peak1.bed -b peak2.bed -wb -filenames \
  | cut -f1,2,3,4 | sort | uniq \
  | cut -f1,2,3 | uniq -c \
  | awk '{sub(/^[ ]+/, "")}' \
  | awk '{print $1}' | sort | uniq -c