當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

Hemberg-lab单细胞转录组数据分析（三）

發布時間：2025/3/15 编程问答 22 豆豆

生活随笔收集整理的這篇文章主要介紹了 Hemberg-lab单细胞转录组数据分析（三）小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

Hemberg-lab單細胞轉錄組數據分析（一）

Hemberg-lab單細胞轉錄組數據分析（二）

收藏|北大生信平臺"單細胞分析、染色質分析"視頻和PPT分享

生信老司機以中心法則為主線講解組學技術的應用和生信分析心得 - 限時免費

scRNA-seq原始數據加工

FastQC

得到單細胞RNA-seq測序數據后，首先檢查測序reads的質量。為了完成這個任務，我們使用的工具是FastQC。FastQC是一款質控工具，能對bulk RNA-seq和單細胞RNA-seq的原始數據進行質量控制 (其他類型的高通量測序結果也適用)。FastQC以原始測序reads為輸入(fastq格式)，輸出序列質量報告。復制粘貼下面的鏈接到你的瀏覽器進入FastQC官網：

https://www.bioinformatics.babraham.ac.uk/projects/fastqc/

這個網址包含下載和安裝FastQC及示例報告文件的鏈接。向下滾動頁面到Example Reports并點擊Good Illumina Data，會看到高質量Illumina測序Reads的理想質控結果。如果使用Docker鏡像，則FastQC已經安裝好。如果是自己的服務器，FastQC下載下來即可使用（依賴于Java）。生信寶典的推文NGS基礎 - FASTQ格式解釋和質量評估對FASTQ原始數據和FastQC的使用和結果描述有比較詳細的介紹，如果不熟悉，建議閱讀。

現在讓我們自己來做一個FastQC報告。

今天我們會用一個mESC數據集 [@kim_characterizing_2015] (原文檔文獻引用錯誤，有更新)的單細胞測序結果為例做后續分析。細胞使用SMART-seq2方法構建測序文庫并進行雙端測序。這些文件在Share文件夾下。

注意: 當前課程教案是寫給線下參加我們課程的學員，程序和數據都已放置在特定位置。如果您只是跟著教程學習，需要自己配置服務器和原始數據文件，比如下載文件（ERR522959_1.fastq?和ERR522959_2.fastq)并放置到Share?(注意大小寫)目錄。在EBI的ArrayExpress搜索E-MTAB-2600獲取數據下載鏈接：https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-2600/。

NCBI SRA：?NGS基礎 - 測序原始數據下載查看如何從NCBI下載原始測序數據。

現在讓我們來看看文件 (基本Linux命令見免費Linux系統和生信寶典原創學習教程)：

# 很常用的查看文件的命令 # 按 q 退出，空格翻頁 # 更多見上面的鏈接 less Share/ERR522959_1.fastq less Share/ERR522959_2.fastq

任務1：嘗試找到可以生成FastQC報告的命令。提示：嘗試執行

# program -h 或 program --help是最常用的獲取幫助信息的命令 (通常一個短線后面跟單個字母稱為短參數，兩個短線后面跟數個單詞稱為長參數) # program -help 或 program 或program -?也是偶爾會用到的獲取幫助信息的命令，不常用但有，比如blastn -help, samtools -? fastqc -h

這個命令會告訴你FastQC有哪些選項可以設置。如果您看明白了輸出，也就應該知道如何用FastQC進行質控了。運行成功后，正向和反向reads都會獲得一個?.zip和.html。如果成功了，跳到下一部分吧。

解決方案和下載報告

如果還沒成功，可以用下面的命令生成FastQC報告：

# 一般加-p，文件夾存在也不會報錯 # man mkdir查看下-p是什么意思 mkdir -p fastqc_results fastqc -o fastqc_results Share/ERR522959_1.fastq Share/ERR522959_2.fastq

命令執行完成后，可以得到四個文件—每端reads對應一個zip文件和一個html文件。報告在html文件里，我們可以用filezilla或者scp把它從服務器下載到你的電腦上來查看它。如果使用的是Rstudio，在右下角窗口直接打開查看即可。

大概瀏覽一下文件，記住要看雙端reads的報告！什么樣的reads質量好呢？有什么我們需要在意的問題呢？我們該怎么解決這些問題？更多請閱讀生信寶典的推文NGS基礎 - FASTQ格式解釋和質量評估一文。

如果文件很多，一個個看比較麻煩，可以試試整合QC質控結果的利器——MultiQC。

移除接頭和低質量堿基

TrimGalore!工具是對cutadapt軟件的封裝，用于移除測序接頭序列或測序序列末端低質量堿基。

鑒于上一步FastQC報告中有一些接頭污染，因此需要去除接頭部分。看下NGS基礎 - 高通量測序原理，想想什么時候會有接頭產生？

任務2：我們的測試數據中出現了哪種類型接頭？提示：看一下FastQC報告中的Adapter Content部分。

現在讓我們嘗試用Trim Galore!去除那些麻煩的接頭序列，去除接頭后，最好再運行FastQC查看接頭是否去除干凈。

任務3：嘗試自己寫出去接頭的命令提示1：你可以用

查看下Trim Galore的參數描述。

提示2：仔細閱讀上面命令的執行結果 (在最開始使用每個軟件前，都需要仔細閱讀其幫助文檔)，鑒于序列中出現的接頭十分常見，你覺得去除接頭時需要知道接頭具體序列嗎？

任務3：對trimmed reads進行再次FastQC質控分析，接頭序列都去除干凈了嗎？

一旦你覺得已經成功去除接頭序列并用FastQC評估確認了，你可以在下一節核對你的結果。

解決方案

你可以用下面的命令來去除Nextera序列接頭：

# 注意路徑問題 mkdir -p fastqc_trimmed_results trim_galore --nextera -o fastqc_trimmed_results Share/ERR522959_1.fastq Share/ERR522959_2.fastq

記住為trimmed reads文件生成新的FastQC報告！在報告里應該可以看到你的reads的Adaptor Content結果為PASS了。

祝賀！現在已經生成了reads質量報告并實行了接頭去除，下一步，我們會用STAR(STAR有soft-clip機制，理論上只要文庫質量不太差，不進行質控也可以)和Kallisto把質控后的reads比對到參考轉錄組。

文件格式

生信寶典的生信分析過程中這些常見文件的格式以及查看方式你都知道嗎？有比較詳細的介紹，建議跟下面的內容結合著看。

FastQ

FastQ是scRNASeq數據中最原始數據的格式。所有的scRNA-seq方案都會進行雙端測序，根據文庫構建方法不同，條形碼序列 (barcodes)可能出現在測序的左端或右端序列。但是使用唯一分子標簽 (UMIs)的測序方案會產生包含細胞和UMI barcode再加接頭序列但沒有轉錄本序列的reads。因此reads雖然是雙端測序，但比對時按單端reads對待。

FastQ文件格式如下 (NGS基礎 - FASTQ格式解釋和質量評估對FASTQ原始數據有比較詳細的介紹，如果不熟悉，建議閱讀):

>ReadID READ SEQUENCE + SEQUENCING QUALITY SCORES

BAM

BAM文件是存儲比對結果的標準有效的高度壓縮的二進制格式，其文本格式SAM是直接可讀的。BAM/SAM文件包含一個頭部，記錄樣本準備、測序和比對的信息；后面是每個reads的比對結果，tab作為列分隔符。

比對行有下列標準格式：

QNAME：read編號（如果是UMI文庫，通常包含UMI條形碼）

FLAG：數字標記指示reads比對的”類型”，如reads是否比對上，是否為properly paired等，Picard網站可以在”類型”和對應的數字之間進行轉換，有更詳細闡述。

RNAME：參考序列編號（比如比對到的染色體名字）。

POS：最左邊的比對位置。

MAPQ：比對質量

CIGAR：表示reads的匹配/不匹配部分 (可能包括soft-clipping)

RNEXT：mate/next reads比對到的參考序列編號

PNEXT：mate/next reads比對到的第一個堿基位置

TLEN：模板長度（read比對到的參考區域的長度）

SEQ：read序列信息

QUAL：read質量信息

BAM/SAM文件可以用samtools互相轉換。

# 新版本samtools中-S選項忽略，不需要再加，會自己判斷輸入的是bam還是sam格式 samtools view -S -b file.sam > file.bam # -h:包含header samtools view -h file.bam > file.sam

一些測序服務機構會自動把測序reads比對到標準基因組并提供BAM或CRAM格式文件，通常這些基因組不會包含ERCC序列，繼而不會有ERCC reads比對到BAM/CRAM文件中。為了量化ERCCs（或序列有任何其他遺傳變異或外源表達需要考慮時），或者如果想使用不同于標準流程（通常是過時的流程，所以要參加易生信單細胞轉錄組專題培訓）中的比對算法，那你需要把BAM/CRAM文件轉回FastQ:

bedtools可以把BAM文件轉成FastQ，為了避免把比對到多個基因組位置的同一個reads轉換為FASTQ中的多條reads，需要先把BAM文件按read名稱排序，并使用Samtools刪除次級比對 (secondary alignments)。此外Picard里也有把BAM轉成FastQ文件的方法。

# sort reads by name (-n) samtools sort -n original.bam -o sorted_by_name.bam # remove secondary alignments (-F 256) samtools view -b -F 256 sorted_by_name.bam -o primary_alignment_only.bam # convert to fastq bedtools bamtofastq -i primary_alignment_only.bam -fq read1.fq -fq2 read2.fq

bedtools是高通量測序分析的常用工具集，Bedtools使用簡介有一些介紹，可以了解下。

CRAM

CRAM文件與BAM文件類似，只是header中包含比對用到的參考基因組的信息。這使得每個read中與參考基因組相同的堿基可以被進一步壓縮。CRAM也支持一些有損數據壓縮方式來進一步優化儲存，CRAMs格式主要是Sanger/EBI測序機構在使用。

CRAM和BAM文件可以用最新一版的samtools(>=v1.0)相互轉換。然而這個轉換需要預先下載并緩存參考基因組。可以從CRAM文件的頭部元數據中獲得參考基因組信息自行下載，具體按下面的操作進行轉換：

# 設置環境變量，指定基因組文件放置的位置 export REF_CACHE=/path_to/cache_directory_for_reference_genome # -T后面直接跟基因組文件名字；如果上一步沒設置，直接給全路徑也可以 samtools view -b -h -T reference_genome.fasta file.cram -o file.bam samtools view -C -h -T reference_genome.fasta file.bam -o file.cram

手動查看文件

一些時候，需要自己查看文件 (初學時，一定多看看這些文件格式，理解每一部分的含義和參數的使用)，比如查看下文件的header信息是否正確。less和more可以用來在命令行查看任意大小文本文件 (個人更常用less)。管道符|可以在多個命令之間傳輸數據，省卻把中間數據存儲多個拷貝的過程，既簡潔又快速。具體見生信寶典出品的Linux - 管道、標準輸入輸出教程。

less file.txt more file.txt # counts the number of lines in file.txt wc -l file.txt samtools view -h file.[cram/bam] | more # counts the number of lines in the samtools output samtools view -h file.[cram/bam] | wc -l

練習

假如你已經有了一個小的cram格式文件：EXAMPLE.cram

任務1：這個文件是怎么獲得的？用了什么比對軟件?基因組的版本是什么？（提示：檢查頭部文件）
任務2: 有多少reads比對或沒比對上？文件共有多少短序列？secondary alignments有多少？（提示：用FLAG）
任務3：將CRAM格式轉成Fastq文件。轉換后的每條read都是只有一個拷貝嗎？（轉換后的文件命名為?10cells_read1.fastq?10cells_read2.fastq）

小技巧：如果運行某個軟件的幫助命令卡住時，直接輸入命令，看看有無提示信息？

答案

samtools view -T 2000_reference.transcripts.fa -H EXAMPLE.cram | less samtools view -T 2000_reference.transcripts.fa -f 4 EXAMPLE.cram | wc -l # unmapped samtools view -T 2000_reference.transcripts.fa -F 260 EXAMPLE.cram | wc -l # mapped samtools view -T 2000_reference.transcripts.fa -F 256 EXAMPLE.cram | wc -l # total samtools view -T 2000_reference.transcripts.fa -f 256 EXAMPLE.cram | wc -l # secondary alignmentssamtools view -b -h -T 2000_reference.transcripts.fa EXAMPLE.cram -o EXAMPLE.bam samtools sort -n EXAMPLE.bam -o sorted_EXAMPLE.bam samtools view -b -F 256 sorted_EXAMPLE.bam -o primary_EXAMPLE.bam # convert to fastq bedtools bamtofastq -i primary_EXAMPLE.bam -fq 10cells_read1.fq -fq2 10cells_read2.fq

基因組和基因注釋（FASTA,GTF）

為了進行序列比對，需要參考基因組和/或基因組注釋文件（GTF格式或GFF格式）。模式生物的基因組可以從任何主要的數據庫下載：Ensembl, NCBI,或者UCSC Genome Browser。具體操作見NGS基礎 - 參考基因組和基因注釋文件。

GTF文件有基因、轉錄本和外顯子的注釋 (NGS基礎 - GTF/GFF文件格式解讀和轉換中描述更詳細具體)，一個9列的文件：

seq_id：序列的編號，一般為chr或scafold編號；

source: 注釋的來源，一般為數據庫或者注釋的機構，如果未知，則用點.代替；

type: 注釋信息的類型，比如Gene、cDNA、mRNA、CDS等;

start: 該基因或轉錄本在參考序列上的起始位置 (從1開始，包含);

end: 該基因或轉錄本在參考序列上的終止位置 (從1開始，包含);

score: 得分，數字，是注釋信息可能性的說明，可以是序列相似性比對時的E-values值或者基因預測是的P-values值，.表示為空;

strand: 該基因或轉錄本位于參考序列的正鏈(+)或負鏈(-)上;

phase: 僅對注釋類型為CDS有效，表示起始編碼的位置，有效值為0,1,2。(對于編碼蛋白質的CDS來說，本列指定下一個密碼子開始的位置。每3個核苷酸翻譯一個氨基酸，從0開始，CDS的起始位置，除以3，余數就是這個值，表示到達下一個密碼子需要跳過的堿基個數。0表示該編碼框的第一個密碼子第一個堿基位于其5’末端；1表示該編碼框的第一個密碼子的第一個堿基位于該編碼區外；2表示該編碼框的第一個密碼子的第一、二個堿基位于該編碼區外；)；

attributes: 一個包含眾多屬性的列表，格式為標簽值(tag value)，以多個鍵值對組成的注釋信息描述，鍵與值之間用`` (空格分割)，不同的鍵值用;隔開，如gene_id “gene”; transcript_id “geneA.1”; database_id “0012”; modified_by “Damian”; duplicate 0。

根據我們的經驗，Ensembl是最容易使用的，并且具有最大的注釋集。NCBI往往更嚴格，僅包含高可信度的基因注釋。而UCSC包含多個使用不同標準的基因注釋。

如果你的實驗系統含有非標準序列 (比如ERCC spike-ins)，那么這些序列必須加到基因組fasta和gtf上來定量它們的表達。另外，對CRISPR相關序列或其他過表達/報告載體也必須進行相同的操作。

為了獲得最大的可用性和靈活性，我們建議為任何添加的非標準序列創建完整詳細的fasta序列文件和gtf文件。

目前還沒有標準化的方法來做到這一點。下面是我們寫的一個perl腳本 (和生信寶典寫的awk腳本)，可以把ERCC序列轉成對應的gtf和fasta文件，以便附加到標準基因組中。如果要量化內含子reads時，您可能還需要更改gtf文件以處理內含子中的重復元件。任何腳本語言甚至是awk或一些文本編輯器都可以相對有效地完成這項任務，但它們超出了本次課程的范圍。(具體教程之前都有過推文介紹：awk學習見Linux - 常用和不太常用的實用awk命令，python學習見Python極簡教程（一）。對應的視頻見http://bioinfo.ke.qq.com/。)

# 下面兩個awk命令可以實現perl腳本的功能 # 如果是windows的文件，替換下末尾的換行符 sed -i 's/^M//' ERCC_Controls_Annotation.txt # 轉成FASTA序列 awk 'BEGIN{OFS=FS="\t"}{if(FNR>1) print ">"$1"\n"$NF"NNNN"}' ERCC_Controls_Annotation.txt >ERCC.fa # 轉成GTF awk 'BEGIN{OFS=FS="\t"}{if(FNR>1) {seq_len=length($NF)+2; attr="gene_id \""$1"-"$2"\"; transcript_id \""$1"-"$2"\"; exon_number \"1\"; gene_name \""$1"-"$2"\""; print $1,"ERCC","gene",1,seq_len,".","+",".",attr; print $1,"ERCC","transcript",1,seq_len,".","+",".",attr; print $1,"ERCC","exon",1,seq_len,".","+",".",attr;}}' ERCC_Controls_Annotation.txt >ERCC.gtf

下面的perl腳本存儲為文件?erccToAnno.pl，并與ERCC_Controls_Annotation.txt置于同一目錄下，運行perl erccToAnno.pl會生成兩個文件ERCC_Controls.fa和ERCC_Controls.gtf，即為結果。與上面生成的ERCC.fa和ERCC.gtf一樣。(awk腳本中基因名字做了修改，不允許空格出現)

# Converts the Annotation file from # https://www.thermofisher.com/order/catalog/product/4456740 to # gtf and fasta files that can be added to existing genome fasta & gtf files. # 下面的perl腳本存儲為文件 erccToAnno.pl，并與ERCC_Controls_Annotation.txt置于同一目錄下， # 運行perl erccToAnno.pl會生成兩個文件ERCC_Controls.fa和ERCC_Controls.gtf。my @FASTAlines = (); my @GTFlines = (); open (my $ifh, "ERCC_Controls_Annotation.txt") or die $!; <$ifh>; #header while (<$ifh>) {# Do all the important stuffchomp;my @record = split(/\t/);my $sequence = $record[4];$sequence =~ s/\s+//g; # get rid of any preceeding/tailing white space$sequence = $sequence."NNNN";my $name = $record[0];my $genbank = $record[1];push(@FASTAlines, ">$name\n$sequence\n"); # is GTF 1 indexed or 0 indexed? -> it is 1 indexed # + or - strand?push(@GTFlines, "$name\tERCC\tgene\t1\t".(length($sequence)-2)."\t.\t+\t.\tgene_id \"$name-$genbank\"; transcript_id \"$name-$genbank\"; exon_number \"1\"; gene_name \"ERCC $name-$genbank\"\n");push(@GTFlines, "$name\tERCC\ttranscript\t1\t".(length($sequence)-2)."\t.\t+\t.\tgene_id \"$name-$genbank\"; transcript_id \"$name-$genbank\"; exon_number \"1\"; gene_name \"ERCC $name-$genbank\"\n");push(@GTFlines, "$name\tERCC\texon\t1\t".(length($sequence)-2)."\t.\t+\t.\tgene_id \"$name-$genbank\"; transcript_id \"$name-$genbank\"; exon_number \"1\"; gene_name \"ERCC $name-$genbank\"\n"); } close($ifh);# Write output open(my $ofh, ">", "ERCC_Controls.fa") or die $!; foreach my $line (@FASTAlines) {print $ofh $line; } close ($ofh);open($ofh, ">", "ERCC_Controls.gtf") or die $!; foreach my $line (@GTFlines) {print $ofh $line; } close ($ofh);

總結

以上是生活随笔為你收集整理的Hemberg-lab单细胞转录组数据分析（三）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。