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入れないとリード名がオリジナルと違うものになってしまう (赤枠部分が追加されている)
![sra.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/383410/935b63a2-7c15-8d77-ea25-fb303311385c.png)

本来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を除く。(本当は残しておきたいがとりあえず消す)

```bash

$ sed -e 's/barcode=[^-]*-[^-]*\-[^-]*\-[^-]*\-/barcode=/' test.fq > modified.fq
```
![modified.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/383410/fbe1af0d-bd5d-647f-bd9f-2fefadf962ea.png)

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

```bash

$ 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
```
###インデックス作成
```bash

$ 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
```
#マッピング
```bash

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

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

#SAMファイル処理
SAMファイルでは12行目以降にタグとして追加情報を記載できる。今回はマッピング結果のSAMファイルにバーコードとUMI情報のタグを追加する。[この形式](https://github.com/samtools/hts-specs/blob/master/SAMv1.pdf)に従います。

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

```bash

$ 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を合わせて入れておく。

```r: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に変換し、インデックスを作成

```bash
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, シーケンスエラーを加味してくれる。

```bash
$ 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 #インプットとアウトプット
```
#バーコードをウェルに変換する
[ここ](ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE107nnn/GSE107947/suppl/GSE107947_experimental_design.txt.gz)でPB, CBとウェルIDの対応表がダウンロードできるので、これを元にウェルIDに変換する

```R
#ウェルの表とカウントデータ読み込み
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で変換

```R
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)

```

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