劉小澤寫于19.6.13-第二單元第一講：單細胞轉(zhuǎn)錄組上游分析之shell回顧

筆記目的：根據(jù)生信技能樹的單細胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

前言

目前主流的單細胞測序技術(shù)主要有兩種：主打基因數(shù)量的smart-seq2和主打細胞數(shù)量的10X Genomics。單細胞轉(zhuǎn)錄組分析和常規(guī)的轉(zhuǎn)錄組分析沒有太大區(qū)別，只是將原來作為一個樣本的一塊組織給分解，變成大量細胞，并且每個細胞單獨作為一個樣本；就像TCGA這樣的大型人群隊列中測一千個人的轉(zhuǎn)錄組，只不過一次性將這一千個人的轉(zhuǎn)錄組放在一起進行分析。

單細胞數(shù)據(jù)的一個特點就是：每個樣本的數(shù)據(jù)量小。以人為例，常規(guī)轉(zhuǎn)錄組一般都能測30M、50M，也就是動輒幾千萬條reads。但是單細胞能夠測500萬條reads就非常厲害了

單細胞上游分析需要get的點

需要用到linux、R以及常規(guī)轉(zhuǎn)錄組分析流程

拿數(shù)據(jù)

首先拿到文章，先搜索"GEO"或者"GSE"

然后點進去超鏈接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111229，就到了GEO數(shù)據(jù)庫，注意三點就可以：

• 只需要改鏈接最后的GSExxxxx就可以快速訪問不同的GEO數(shù)據(jù)

• 要是想下載作者做好的表達矩陣，然后直接進行下游分析

  
  

• 要是想從原始測序數(shù)據(jù)fq文件開始自己分析，就要進入SRA測序數(shù)據(jù)中心

  
  

進入測序數(shù)據(jù)中心就可以看到：接近800個細胞才測了不到10G的數(shù)據(jù)量，要知道一個常規(guī)轉(zhuǎn)錄組數(shù)據(jù)都可以做到10G數(shù)據(jù)量 (另外，10X的數(shù)據(jù)量要比smart-seq2還要小)。

下載+轉(zhuǎn)換

利用conda安裝sra-tools 就會帶prefetch ，然后利用這個軟件下載sra原始數(shù)據(jù)，默認使用http方式下載，如果網(wǎng)速比較慢，還可以使用"開掛"模式，具體操作看：來吧，加速你的下載

下載全部的sra文件后會發(fā)現(xiàn)：smart-seq2結(jié)果中每個細胞都是一個單獨的sra文件，它是單細胞的單樣本

這一點和10X是有區(qū)別的，10X是一個樣本中就包含了4000-8000個細胞，但不會拆分成4000-8000個fq文件，需要進一步利用UMI、barcode將細胞分開，也就是說，10X多了一步拆分的過程。

有了sra，就需要轉(zhuǎn)換成fastq，利用的軟件就是fastq-dump ，它也屬于sra-tools ，一般參數(shù)設定為：fastq-dump --gzip --split-3 -O ./

因為采用了3’單端測序，因此轉(zhuǎn)換后也只有一個fq文件(如果是采用10X，雖然也是單端測序，但是它轉(zhuǎn)換結(jié)果會是三個文件：sample index、barcode+UMI、真正的測序reads)

拿到原始數(shù)據(jù)后一般要利用fastqc和multiqc進行質(zhì)控，質(zhì)量合格進入下一環(huán)節(jié)

比對

利用hisat2 ，需要注意：先利用hisat2-build命令構(gòu)建基因組索引，文章使用的是小鼠的基因組mm10

只要有一臺服務器存在這個索引即可，可以使用跨服務器拷貝scp命令
使用方法簡單：
https://linuxtools-rst.readthedocs.io/zh_CN/latest/tool/scp.html

scp local_file remote_username@remote_ip:remote_folder
# 拷貝整個目錄就加參數(shù)-r

# 因為是單端測序文件，因此hisat2軟件就要用-U選項(如果雙端的話，直接用 -1 -2 選項即可)
index=/PATH_to_hisat2_mm10/
ls *.gz | while read i;do hisat2 -p 10 -x $index -U $i -S ${i%%.*}.hisat.sam;done
# -S選項指的是輸入比對結(jié)果SAM文件，它的參數(shù)就是輸出的SAM叫什么?？吹狡渲杏幸粋€%%.*，它的意思就是取我們輸入原始測序文件(i)的前綴

直接運行的話，會給出比對結(jié)果，看一下比對率

得到的結(jié)果：

比對后SAM要轉(zhuǎn)為BAM才能進行下一步定量，利用samtools

ls *.sam|while read i ;do (samtools sort -O bam -@ 10  -o $(basename ${i} ".sam").bam ${i});done
# 其中一個小tip就是：basename這里，會返回每個$i，也就是每個sam文件的名稱，然后".sam"就是將名稱中的.sam去掉，于是只留下了前綴名，這樣才能更改文件名為bam

可以看到bam比sam小了很多倍

轉(zhuǎn)換后繼續(xù)構(gòu)建索引：

ls *.bam |xargs -i samtools index {}

定量

它的目的就是將bam比對結(jié)果和參考的基因注釋(其實也是基因在基因組上的位置信息)進行比較，看看我們的結(jié)果中對應了哪些基因。主要利用featureCounts

GTF文件可以去Gencode數(shù)據(jù)庫下載：
ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M21/gencode.vM21.annotation.gtf.gz

# 還是先指定gtf的位置
gtf=/YOUR_GTF_PATH/
featureCounts -T 5 -t exon -g gene_id  -a $gtf -o all.id.txt  *.bam  1>counts.id.log 2>&1 &
# 如果是雙端測序，可以加上-p選項

關于這款軟件：featureCounts官網(wǎng)

• featureCounts不僅支持基因或轉(zhuǎn)錄本的定量，還可以進行exons, promoter, gene bodies, genomic bins and chromosomal locations的定量

• 這個軟件不是單獨下載的，它是集成在subread軟件中，因此只要下載好subread就能使用featureCounts。

• 它需要的輸入文件也很簡單：比對的sam/bam文件(我們經(jīng)常使用bam，是因為它占據(jù)硬盤空間小)、注釋文件GTF

• 它的定量有兩個層次：一個是對feature定量，另一個是對metafeature 進行定量。官網(wǎng)對它們定義的描述：

  Each entry in the provided annotation file is taken as a feature (e.g. an exon). A meta-feature is the aggregation of a set of features (e.g. a gene).

  feature也就是基因組區(qū)間最小的信息(如外顯子)；metafeature可以是多個feature的組合，如同一個基因的多個外顯子集合；因此這款軟件可以單獨對外顯子定量，也可以對基因進行定量。但只有比對到多個不同區(qū)間時，才會分別計數(shù)

• -T表示線程數(shù)，默認是1；-t 表示要計數(shù)的feature名稱，也就是GTF的第三列信息，默認是exon；-g表示提取的GTF最后一列attribute信息，默認是gene_id；-a 指定使用的注釋文件；-o是輸出文件

得到的結(jié)果中可以看到，基因?qū)耐怙@子、起始終止坐標等。其中大量的基因在7個示例樣本中都沒有表達量，也就是說，對于大部分細胞來說許多基因都是測不到的

得到的表達矩陣也就是文中一開始提到的作者給出的rawcounts表達矩陣，只不過這里我們只有7個樣本，而真正有768個樣本。另外表達矩陣的軟件、版本、參數(shù)有所差別，因此得到的不會完全一樣，這也是可接受的