scRNA-seqデータの再解析記録 (MARS-seq)(未)

fastqから再解析する

使うデータ

GSM2884123 (https://www.ncbi.nlm.nih.gov/sra/?term=SRR6366715)
マウスにH1N1感染させて肺のCD45+細胞をMARS-seqで解析したもの

MARS-seq

MARS-seqについては[元論文のSupplement] (https://science.sciencemag.org/highwire/filestream/595372/field_highwire_adjunct_files/0/Jaitin-SM.pdf) か[このサイト] (https://teichlab.github.io/scg_lib_structs/methods_html/MARS-seq.html)が参考になった。

今後はPool BarcodeをPB, Cell BarcodeをCBと書きます。

fastqダウンロード

2ファイルに分かれているのでfastq-dumpしてマージする。

$ fastq-dump --split-files --origfmt --defline-qual '+' SRR6366715
$ fastq-dump --split-files --origfmt --defline-qual '+' SRR6366716
$ cat SRR6366715_1.fastq SRR6366716_1.fastq > test.fastq 

※SRAからダウンロード時には、--origfmt入れないとリード名がオリジナルと違うものになってしまう (赤枠部分が追加されている)

本来MARS-seqではトランスクリプトのReadとCB,UMIのReadに分かれているが、今回のfastqはリード名の後にバーコードのデータが付加されている。
(barcode=tail-PB_Phred-CB_Phred-UMI_Phred-PB-CB-UMI)

他の方式のscRNA-seqでも同様になっているものがあったが、何のアプリケーションでバーコード抽出したらこうなるのか、よく分からなかった。(bcl2fastqならコロン区切り、UMI-toolsならアンダースコア区切りで出力されるはず...)

前処理する

STARでマッピングする時にリード名が/で切られてしまうので、リード名からPhredを除く。(本当は残しておきたいがとりあえず消す)

$ sed -e 's/barcode=[^-]*-[^-]*\-[^-]*\-[^-]*\-/barcode=/' test.fq > modified.fq

マッピング準備

リファレンスをダウンロード

GRCm38のmm10からゲノムFastaとChromosomeごとの遺伝子GTFをダウンロード。
NCBIとEnsemblの2種類あるが、NCBIの方には最後の解析で使うミトコンドリアDNAのアノテーションデータが入ってないので、Ensemblからダウンロードする。

$ wget ftp://ftp.ensembl.org/pub/release-99/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
$ gunzip Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
$ wget ftp://ftp.ensembl.org/pub/release-99/gtf/mus_musculus/Mus_musculus.GRCm38.99.chr.gtf.gz
$ gunzip Mus_musculus.GRCm38.99.chr.gtf.gz

インデックス作成

$ mkdir STAR_index_mouse
$ STAR \
--runThreadN 12 \
--runMode genomeGenerate \
--genomeDir ./STAR_index_mouse \ #インデックス作るディレクトリ指定
--genomeFastaFiles ./Mus_musculus.GRCm38.dna.primary_assembly.fa \ #ゲノムFasta
--sjdbGTFfile ./Mus_musculus.GRCm38.99.chr.gtf #GTF

マッピング

$ STAR \
--runThreadN 12 \
--genomeDir ./STAR_index_mouse \ #リファレンス
--readFilesIn modified.fq \ #データ入力
--outReadsUnmapped Fastx \ #マッピングされなかったリードをFastqで出力
--outSAMunmapped Within \ #マッピングされなかったリードもSAMに含める
--quantMode TranscriptomeSAM GeneCounts #転写産物ベースのマッピング結果と遺伝子カウントも出力

今回は転写産物ベースのマッピング結果を使う。

SAMファイル処理

SAMファイルでは12行目以降にタグとして追加情報を記載できる。今回はマッピング結果のSAMファイルにバーコードとUMI情報のタグを追加する。この形式に従います。

先に@から始まるヘッダーと本体を分けておく。

$ grep "^@" Aligned.toTranscriptome.out.sam > header.txt
$ grep -v "^@" Aligned.toTranscriptome.out.sam > no_header.txt

Tag	中身
BC	PB
CB	CB
RX	UMI
XX	PB + CB

この後使うUMI-toolsがPBにまで対応していないので、XXタグにPBとCBを合わせて入れておく。

Tag_attaching.R

library(stringr)

df <- read.table("no_header.txt",header=FALSE,fill=TRUE,row.names = NULL,check.names=FALSE,comment.char="")
df[,ncol(df)+1]  <- str_sub(df[,1],-21,-1)
df <- filter(df, str_detect(df[,ncol(df)], ".*N.*")==FALSE)
df <- df[!duplicated(df[,ncol(df)]),]
bc <- str_sub(df[,ncol(df)],1,4)
cb <- str_sub(df[,ncol(df)],6,12)
umi <- str_sub(df[,ncol(df)],14,21)

df <- df[,1:11]
df[,ncol(df)+1] <- paste0("BC:Z:",bc)
df[,ncol(df)+1] <- paste0("CB:Z:",cb)
df[,ncol(df)+1] <- paste0("RX:Z:",umi)
df[,ncol(df)+1] <- paste0("XX:Z:",bc,cb)
write.table(df,"Aligned_processed.txt",quote=F,col.names = F,row.names = F,sep="\t")

headerと繋いで、ソートしつつBAMに変換し、インデックスを作成

cat header.txt Tag_attached.txt > Tag_attached.sam
samtools sort -O bam -o Tag_sorted.bam Tag_attached.sam
samtools index Tag_sorted.bam

UMI処理

UMI-toolsを使ってUMIの重複を削除し、細胞単位でリードをカウントする。(UMI-toolsはbiocondaからインストール可能)
UMI-toolsはUMIのPCR, シーケンスエラーを加味してくれる。

$ umi_tools dedup \ #UMIで重複削除
--extract-umi-method=tag --umi-tag=RX \ #UMIタグ(RX)を指定
--per-cell --cell-tag=XX \ #セルバーコード(XX)を指定
--per-gene --per-contig \
--output-stats dedup \ #UMIの統計情報を出力。dedupはファイル名のprefix。
-I Tag_sorted.bam --out-sam -S dedup.sam #インプットとアウトプット

$ umi_tools count \ #細胞ごとにカウント
--extract-umi-method=tag --umi-tag=RX \ #UMIタグ(RX)を指定
--per-cell --cell-tag=XX \#セルバーコード(XX)を指定
--per-gene --per-contig \
-I Tag_sorted.bam -S umi_count.txt #インプットとアウトプット

バーコードをウェルに変換する

ここでPB, CBとウェルIDの対応表がダウンロードできるので、これを元にウェルIDに変換する

# ウェルの表とカウントデータ読み込み
mapped <- read.table("umi_count.txt",fill=T,sep="\t",header=T)
cell.data <- read.table("GSE107947_experimental_design.txt",fill=T,sep="\t",skipNul=T,header=T)
cell.data2 <- data.frame(cell.data$Well_ID,cell.data$Batch_description,paste0(cell.data$Pool_barcode,cell.data$Cell_barcode))
colnames(cell.data2) <- c("Well_ID","Batch_description","barcode")
cell.data2 <- cell.data2[cell.data2$Batch_description == "Influenza_treated_CD45+_48h_replicate_II"]

# カウントとウェルの表をマージ
merged <- merge(cell.data2,mapped,by.x="barcode",by.y="cell",sort=F,all=T)
merged <- na.omit(merged)

# クロス集計
crs <- xtabs(count~gene+Well_ID,data=merged)
crs <- as.data.frame.matrix(crs)

ここでマージしてNAを含む行を削除したときにカウントデータからデータが大量になくなってしまった。
PCR，シーケンスエラーでウェルデータにないバーコードが多く含まれてしまっているのが原因だと思うが、その辺りを補正するツールを見つけられなかった。

TranscriptをGeneSymbolに

biomaRtで変換

library(biomaRt)
db <- useMart("ensembl")
hg <- useDataset("mmusculus_gene_ensembl",mart=db)

res <- getBM(attributes = c("ensembl_transcript_id","external_gene_name"),filters="ensembl_transcript_id",value=rownames(crs),mart=hg)

crs <- cbind(res$external_gene_name,crs)
colnames(crs)[1] <- "gene_name"
crs <- aggregate(x=crs[,2:ncol(crs)],by=list(crs$gene_name),sum)

write.table(crs,"dge.tsv",quote=F,sep="\t",row.names = F)

おわってない

バーコードの補正方法ご存じの方いれば教えてください