サルマップ2018のまとめ（ちょっと賢いループ付き）

2018年06月13日（水）

投稿者: 中川真一

事前準備

OSXインストール

Aspera Connect

Homebrewでインストール（ターミナルから）

/usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"
brew install sratoolkit
brew install samtools
brew install bedtools
brew install wget
brew tap brewsci/bio
brew install hisat2

Rでインストール（Rコンソールから）

source("https://bioconductor.org/biocLite.R")
biocLite(pkgs="Rsubread")
biocLite("DESeq2")
biocLite("Biostrings")
biocLite("rtracklayer")
install.packages("ggplot2")

作業フォルダ（例）

/Users/Nakagawa_imac/Desktop/2018_NGS_standard
/Users/Nakagawa_imac/Desktop/2018_NGS_standard/genome

genomeフォルダにダウンロード&インデックスの作成（ターミナルから）

cd /Users/Nakagawa_imac/Desktop/2018_NGS_standard/genome
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M17/gencode.vM17.chr_patch_hapl_scaff.annotation.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M17/GRCm38.p6.genome.fa.gz
wget https://ncrna.jp/blog/item/download/72_e35f6683783f252239a42685635879e8 #（ダウンロード後拡張子を.bedに変える）
hisat2-build GRCm38.p6.genome.fa genome_index

genomeフォルダの中身の確認

gencode.vM17.chr_patch_hapl_scaff.annotation.gtf
genome_index.1.ht2
genome_index.2.ht2
genome_index.3.ht2
genome_index.4.ht2
genome_index.5.ht2
genome_index.6.ht2
genome_index.7.ht2
genome_index.8.ht2
GRCm38.p6.genome.fa
mm10_rmsk_rRNA.bed

マッピングまで

single end sraファイルのダウンロードとマッピング（SRR6946223からSRR6946228までの連番の場合）（ターミナルから）

#!/bin/sh

##example for single end read sra file download & mapping

#set directory
cd /Users/Nakagawa_imac/Desktop/2018_NGS_standard

#set loop
for srr in SRR69462{23..28}
do
  #sra download
  /Users/Nakagawa_imac/Applications/Aspera\ Connect.app/Contents/Resources/ascp -T -k 1 -i /Users/Nakagawa_imac/Applications/Aspera\ Connect.app/Contents/Resources/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/${srr::6}/${srr}/${srr}.sra .
  #sam to bam conversinon
  fastq-dump --gzip ./${srr}.sra
  #hisat2 mapping
  hisat2 -p 4 -x ./genome/genome_index -U ${srr}.fastq.gz -S ${srr}.sam
  #sam to sorted bam
  samtools sort -@ 4 -O bam -o ${srr}.sort.bam ${srr}.sam
  #create index
  samtools index ${srr}.sort.bam
  #remove rRNA reads
  intersectBed -abam ${srr}.sort.bam -b ./genome/mm10_rmsk_rRNA.bed -v > ${srr}.sort_rRNA_rm.bam
  #remove sam file
done

paired end sraファイルのダウンロードとマッピング（SRR1182778からSRR1182781までの連番の場合）（ターミナルから）

#!/bin/sh

##example for paired end read sra file download & mapping

#set directory
cd /Users/Nakagawa_imac/Desktop/2018_NGS_standard

#set loop
for srr in SRR11827{78..81}
do
  #sra download
  /Users/Nakagawa_imac/Applications/Aspera\ Connect.app/Contents/Resources/ascp -T -k 1 -i /Users/Nakagawa_imac/Applications/Aspera\ Connect.app/Contents/Resources/asperaweb_id_dsa.openssh anonftp@ftp.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/${srr::6}/${srr}/${srr}.sra .
  #sam to bam conversinon
  fastq-dump --gzip --split-files ./${srr}.sra
  #hisat2 mapping
  hisat2 -p 4 -x ./genome/genome_index -1 ${srr}_1.fastq.gz -2 ${srr}_2.fastq.gz -S ${srr}.sam
  #sam to sorted bam
  samtools sort -@ 4 -O bam -o ${srr}.sort.bam ${srr}.sam
  #create index
  samtools index ${srr}.sort.bam
  #remove rRNA reads
  intersectBed -abam ${srr}.sort.bam -b ./genome/mm10_rmsk_rRNA.bed -v > ${srr}.sort_rRNA_rm.bam
done

フォルダ内のシングルリードの.gzファイルを一括処理（ターミナルから）

#!/bin/sh
#set directory
cd /Users/Nakagawa_imac/Desktop/2018_NGS_standard

for file in *.gz; do
  id=${file%.fastq.gz}
  #hisat2 mapping
  hisat2 -p 4 -x ./genome/genome_index -U $file -S $id.sam
  #sam to bam conversinon
  samtools sort -@ 4 -O bam -o $id.sort.bam $id.sam
  #create index
  samtools index $id.sort.bam
  #remove rRNA reads
  intersectBed -abam $id.sort.bam -b ./genome/mm10_rmsk_rRNA.bed -v > $id.sort_rRNA_rm.bam
done

bashのfor文でよく使うものメモのページとを参考にさせていただきました。便利！

速攻でコメントいただきました。すごい。どんどん洗練されていく。ありがとうございます。

普通は for gzfile in *.gz; do とすべきでは。ファイル名にスペースなどかあるとおかしくなるし…。 https://t.co/QOmaV1wSeS
— Hoshi Takanori (@hoshi_takanori) June 18, 2018

これ意図がわからないが単に for x in *.gz; do としない理由は何だろう？ for x in $(ls|grep \.gz$) ではファイル名に空白が含まれていると悲惨なことになりそうな気がする https://t.co/vCCi3AlKKd
— Haruhiko Okumura (@h_okumura) June 19, 2018

マッピング以降

RsubreadとDESeqを用いた発現解析（Rコンソールから）

working directyの中の.sort_rRNA_rm.bamがついたファイル名をlist.filesで取得して集めてまとめて処理してます。別グループのファイルがある時はそれらを別のフォルダに移します。

#count data
library(Rsubread)
setwd("~/Desktop/2018_NGS_standard")
fl <- list.files(pattern = "*.sort_rRNA_rm.bam")
fc <- featureCounts (files = fl,
                     annot.ext="./genome/gencode.vM17.chr_patch_hapl_scaff.annotation.gtf",
                     isGTFAnnotationFile=TRUE,
                     nthreads=4,
                     allowMultiOverlap = TRUE)
write.table(file = "counts.txt", fc$counts, quote=FALSE, sep="\t")

#normalize data/fold change analyses
library(DESeq2)
dt <- as.matrix(fc$counts)
dt1 <- dt[apply(dt,1,sum)>6,]
group <- data.frame(con = factor(c("normal","normal","normal","HFD","HFD","HFD")))
dds <- DESeqDataSetFromMatrix(countData = dt1, colData = group, design = ~ con)
dds <- DESeq(dds)
res <- results(dds)
res$gene_id <- row.names(res)
write.table(res,file="DESeq2_result.txt",sep="\t",row.names=F,col.names=T,quote=F)

gene symbol付きのファイルの作成（Rコンソールから）

#count data
library(rtracklayer)
gtf <- readGFF("./genome/gencode.vM17.chr_patch_hapl_scaff.annotation.gtf")
gtf_gene <- subset(gtf, gtf$type == "gene")
gtf_gene_1 <- gtf_gene[,c("gene_id","gene_name","gene_type")]
DEseq2_result <- data.frame(res)
DESeq2_result_ref <- merge(DEseq2_result,gtf_gene_1)
write.table(res,file="DESeq2_result_ref.txt",sep="\t",row.names=F,col.names=T,quote=F)

ggplot2を用いたMA-plotの作成（Rコンソールから）

#draw MA-plot
#ggplot2
library(ggplot2)
#select subset of genes
all <- DESeq2_result_ref
mRNA <- DESeq2_result_ref[DESeq2_result_ref$gene_type == "protein_coding",]
lincRNA <- DESeq2_result_ref[DESeq2_result_ref$gene_type == "lincRNA",]
#drow MA plot
g <- ggplot(NULL)
g <- g + geom_point (data = all, aes (x = log(baseMean,base = 10), y = log2FoldChange), color = "gray")
g <- g + geom_point (data = mRNA, aes (x = log(baseMean,base = 10), y = log2FoldChange), color = "coral", alpha = 0.5)
g <- g + geom_point (data = lincRNA, aes (x = log(baseMean,base = 10), y = log2FoldChange), color = "springgreen", alpha = 0.5)
#modify plot
g <- g + geom_hline(yintercept = 0,colour="red",alpha = 0.5)
g <- g + theme(axis.text.x = element_text(size = 18),axis.title.x=element_text(size=20)) + #X label size
         theme(axis.text.y = element_text(size = 18),axis.title.y=element_text(size=20)) + #Y label size
         theme(plot.title = element_text(size=22)) + #title size
  ylab("log2 (HFD/control) ") + # X label
  xlab("log10 (baseMean)") +# Y label
  labs(title = "gene expression changes upon HFD feeding")  #title
ggsave(filename = "MA_plot.png",g,width = 5, height = 5)

カテゴリ: 技術紹介

中川真一

北海道大学　薬学研究院　教授
▶ プロフィールはこちら

このカテゴリをもっと見る « サルマップ2018 (3) DESeq2による標準化と可視化までフォルダ内のファイル名を取得してまとめて処理 »

トップへ戻る

サルマップ2018のまとめ（ちょっと賢いループ付き）

事前準備

OSXインストール

Homebrewでインストール（ターミナルから）

Rでインストール（Rコンソールから）

作業フォルダ（例）

genomeフォルダにダウンロード&インデックスの作成（ターミナルから）

マッピングまで

single end sraファイルのダウンロードとマッピング（SRR6946223からSRR6946228までの連番の場合）（ターミナルから）

paired end sraファイルのダウンロードとマッピング（SRR1182778からSRR1182781までの連番の場合）（ターミナルから）

フォルダ内のシングルリードの.gzファイルを一括処理（ターミナルから）

マッピング以降

RsubreadとDESeqを用いた発現解析（Rコンソールから）

gene symbol付きのファイルの作成（Rコンソールから）

ggplot2を用いたMA-plotの作成（Rコンソールから）

中川真一

ブログ記事カテゴリ

ブログアーカイブ

ログイン

サイト内検索

サルマップ2018のまとめ（ちょっと賢いループ付き）

事前準備

OSXインストール

Homebrewでインストール（ターミナルから）

Rでインストール（Rコンソールから）

作業フォルダ（例）

genomeフォルダにダウンロード&インデックスの作成（ターミナルから）

マッピングまで

single end sraファイルのダウンロードとマッピング（SRR6946223からSRR6946228までの連番の場合）（ターミナルから）

paired end sraファイルのダウンロードとマッピング（SRR1182778からSRR1182781までの連番の場合）（ターミナルから）

フォルダ内のシングルリードの.gzファイルを一括処理（ターミナルから）

マッピング以降

RsubreadとDESeqを用いた発現解析（Rコンソールから）

gene symbol付きのファイルの作成（Rコンソールから）

ggplot2を用いたMA-plotの作成（Rコンソールから）

中川 真一

ブログ記事カテゴリ

ブログアーカイブ

ログイン

サイト内検索

中川真一