RNA-seqその５、Hisat2でマッピング

RNAseqの結果をリファレンス配列にマッピングする
mecobalamin.hatenablog.com

レファレンスデータはヒトのゲノムデータ
Human Genome Resources at NCBI - NCBI
ここからReference Genome Sequenceをダウンロード

28 Feb. 2019追記
ヒトゲノムデータにはデータの書き方の違いで2種類あるらしい
GRCh38とhg38の違い(含むミトコンドリア)
chromosomeの記述の仕方に違いがある
JBrowseで表示したり解析に使うにはhg38がいいのかも
Index of /goldenPath/hg38/bigZips

追記ここまで

20 April 2019追記
ヒトゲノムのデータを作り直した
mecobalamin.hatenablog.com
問題あるかもだけどchromosoneのファイルだけ一つにまとめて
インデックスファイルの作成に使用した
以下のスクリプトではname_dbとpostfixを書き換える

追記ここまで

まずヒトゲノムをhisat2-buildでインデックスファイルに変換する

#!/bin/bash

echo "hisat2 build began"

dir_scripts="$HOME/usr/scripts"
dir_working="$HOME/usr/data"
dir_db="${dir_working}/db"
name_db="GRCh38_latest_genomic"
postfix="fna"

hisat2-build ${dir_db}/${name_db}.${postfix} ${name_db}

input file + 数字 + .ht2という形式で複数のファイルが出力される
GRCh38_latest_genomic.1.ht2
今回は８個だったけどいつも同じファイル数か不明

インデックスファイルができたらいよいよマッピング
hisat2の引数にインデックスファイルのパスと
ファイル名の数字の前までを指定する
マージしたFastqもリード1、2のオプションを付けて読み込ませる
samファイルで出力される

マッピングが終わるとsamtoolsでsamからbamにして
ソートしてindexをつける

stringtieでアノテーションを付ける、
とのことなんだけどよくわかってない
出力されたgtfファイルを見ると
bamに含まれる遺伝子名や場所をまとめたファイルのようだ
transcripだけ読めばカウントできそうだけどどうだろう

#!/bin/bash

echo "hisat2 analysis"

echo "set directories"
dir_scripts="$HOME/usr/scripts"
dir_working="$HOME/usr/data"
dir_db="${dir_working}/db"
name_db="GRCh38_latest_genomic"

echo "open list sra"
dir_list_sra=${dir_scripts}
name_sra="list_sra.txt"
list_filename=$(cat ${dir_list_sra}/${name_sra} | sed 's/\.sra//g')

for i in ${list_filename}
  do
    echo "analyze ${i}"
    dir_out="hisat_${i}"

    mkdir -p ${dir_working}/${dir_out}

    echo "hisat began"
    hisat2 \
      -p 4 \
      -x ${dir_db}/${name_db} \
      --rna-strandness RF \
      --dta \
      -1 ${dir_working}/trim_${i}/ps/ps_${i}_gd_1.fastq \
      -2 ${dir_working}/trim_${i}/ps/ps_${i}_gd_2.fastq \
      -S ${dir_working}/${dir_out}/hisat_${i}.sam \

    echo "samtools began"
    samtools view -bS -@ 4 ${dir_working}/${dir_out}/hisat_${i}.sam |\
    samtools fixmate -@ 4 - ${dir_working}/${dir_out}/fixmate_${i}.bam
    samtools sort ${dir_working}/${dir_out}/fixmate_${i}.bam \
      -o ${dir_working}/${dir_out}/sort_${i}.bam \
      -m 10G \
      -@ 4
    samtools index ${dir_working}/${dir_out}/sort_${i}.bam

    echo "stringtie began"
    stringtie \
      ${dir_working}/${dir_out}/sort_${i}.bam \
      -p 4 \
      -o ${dir_working}/${dir_out}/stringtie_${i}.gtf \
      -l stiringtie_${i}

    echo "analysis finished"

  done