B站：RNA-seq轉(zhuǎn)錄組數(shù)據(jù)分析入門實(shí)戰(zhàn)

1 linux常用命令

touch text.txt  #新建文件
rm -rf /var/log/httpd/access #將會刪除/var/log/httpd/access目錄以及其下所有文件、文件夾
rm -f *html #刪除所有html格式文件
rm -f *zip #刪除所有zip格式文件
tar zxvf #解壓tar.gz文件
tar jxvf samtools-1.11.tar.bz2 #解壓tar.bz2文件

#gzip壓縮后的格式為：*.gz，這種壓縮方式不能保存原文件；
gzip buodo #壓縮
gunzip buodo.gz  #解壓

.bashrc中可以添加一些常用快捷方式，也可以將軟件添加到環(huán)境變量
如何用vi編輯和保存文件

vi ~/.bashrc #打開.bashrc冗恨，然后按a進(jìn)行編輯
alias用來設(shè)置快捷方式
比如：
alias ..='cd ..' 返回上一級
alias ..='cd ../..' 返回上上級
設(shè)置完成后按ESC，然后按:wq保存退出。
. .bashrc 啟動bashrc文件且改，剛才設(shè)置的就可以用了。

2 conda的下載和軟件的安裝

conda的下載

復(fù)制miniconda3的lunux版本鏈接地址板驳，-c是支持?jǐn)帱c(diǎn)下載又跛，就是網(wǎng)斷了再次連接可以繼續(xù)下載。
wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

bash Miniconda3-latest-Linux-x86_64.sh #下載好后笋庄，進(jìn)行解壓安裝

vi ~/.bashrc #然后進(jìn)行查看效扫，是否已經(jīng)將miniconda3的bin文件添加到了環(huán)境變量中去。
source .bashrc #執(zhí)行直砂，將miniconda3添加到環(huán)境變量

導(dǎo)入重要的channel菌仁，因為后期會用到
conda config --add channels r
conda config --add channels conda--forge
conda config --add channels bioconda

image.png

軟件安裝，升級和卸載

conda isntall fastqc #安裝fastqc
fastqc -help #判斷是否安裝成功

conda install samtools=1.4.1 #先安裝一個舊版本
which samtools #查看安裝的位置
conda update samtools #查看最新版本并進(jìn)行安裝
conda remove samtools #卸載軟件

編譯安裝軟件

image.png

先下載一個samtools安裝包静暂，鏈接從官網(wǎng)復(fù)制
wget -c https://github.com/samtools/samtools/releases/download/1.11/samtools-1.11.tar.bz2 
tar jxvf samtools-1.11.tar.bz2 #解壓
cd samtools-1.11/  #cd到當(dāng)前軟件目錄下济丘，可以看到一個configure文件，然后進(jìn)行配置洽蛀。

1. 配置(在root下新建一個samtools摹迷，最后將軟件安裝在這里)
 ./configure --prefix='/root/samtools'
報錯處理：
configure: error: no acceptable C compiler found in $PATH 顯示錯誤主要是沒有C編譯器。
安裝C編譯器 :yum -y install gcc

2. 編譯
make
報錯處理： （問題都是一樣郊供，包沒有安裝峡碉，后續(xù)百度安裝即可）
2.1 CentOS編譯安裝軟件過程中遇到zlib.h: No such file or directory，使用命令：yum install zlib-devel  解決問題驮审。
2.2 curses.h: No such file or directory 錯誤解決辦法yum install ncurses-devel ncurses
2.3 centos編譯時報錯：error: bzlib.h: No such file or directory
解釋：bzlib就是bzip2包 .
首先查看你yum源中的bzip2包鲫寄，命令為：yum search bzip2。
然后安裝你yum源中的bzip2包疯淫，執(zhí)行命令：yum install  bzip2-devel.x86_64

3. 安裝
make install

添加環(huán)境變量
安裝完成后地来，如果不在當(dāng)前路徑運(yùn)行samtools報錯，那就需要添加環(huán)境變量

vi ~/.bashrc

image.png

 . ~/.bashrc #加載環(huán)境變量

3 數(shù)據(jù)預(yù)處理

3.1準(zhǔn)備工作

image.png

image.png

image.png

image.png

gunzip *gz #然后在00ref中對兩個gz文件進(jìn)行解壓

3. 原始數(shù)據(jù)的上傳
原始數(shù)據(jù)通過WINscp軟件進(jìn)行上傳。

image.png

#整個操作就在01raw_data文件下進(jìn)行
md5sum *gz >md5.txt #給所有的gz文件生成md5值蜡秽，儲存到md5.txt中府阀。
cat md5.txt #進(jìn)行查看
md5sum -c md5.txt #比對gz文件中已有的md5值和生成的md5是否匹配，匹配返回OK.

image.png

3.2 質(zhì)量控制

image.png

fastqc sample1_R1.fastq.gz #單一樣本質(zhì)控
fastqc sample*gz #批量質(zhì)控開頭為sample载城，gz結(jié)尾的文件
for i in 'ls *gz';do fastqc $i;done #批量質(zhì)控肌似，遍歷該文件夾下所有g(shù)z結(jié)尾的文件，do fastqc诉瓦，然后done結(jié)束

并行批量質(zhì)控

#xargs命令是將上一個命令的輸出作為下一個命令的輸入
#nonhup &是將運(yùn)行程序掛在后臺川队，這樣就可以把每個都掛在后臺，達(dá)到批量的目的睬澡。
#>fastqc.sh是將這個應(yīng)用保存在shell腳本
ls *gz |xargs -I[] echo 'nohup fastqc []&' >fastqc.sh

cat fastqc.sh #可以看到就是將四個程序分別掛到后臺運(yùn)行
#nohup fastqc sample1_R1.fastq.gz&
#nohup fastqc sample1_R2.fastq.gz&
#nohup fastqc sample2_R1.fastq.gz&
#nohup fastqc sample2_R2.fastq.gz&
bash fastqc.sh #運(yùn)行shell腳本
top -c #查看運(yùn)行進(jìn)程固额，ctrl+推出查看

在運(yùn)行的結(jié)果中，可以看到每個樣本有一個html（保存質(zhì)控信息）和zip（保存的文本內(nèi)容）文件

image.png

多個質(zhì)控結(jié)果的查看
首先運(yùn)用conda安裝multiqc煞聪，multiqc支持多種結(jié)果的整合斗躏，同時展示多個結(jié)果文件，可以自動檢測當(dāng)前文件下可以合并的文件昔脯。

multiqc ./

image.png

然后將multiqc_report.html 下載到本地進(jìn)行查看啄糙。

3.3 質(zhì)量過濾

image.png

image.png

接頭序列的選擇：圖中可以看到是Truseq adapter ，所以就根據(jù)單端或者雙端選擇Truseq3云稚。
另外去接頭的參數(shù)選擇TURE,保留下來的reads較多隧饼。

image.png

image.png

用trimmomatic進(jìn)行質(zhì)量過濾，在存在fastq.gz的文件下進(jìn)行運(yùn)行

#在01raw_data下運(yùn)行
# PE是雙端測序静陈，threads 是線程燕雁，../02clean_data/是生成文件的保存位置
#sample1_paired_clean是生成的成對read的文件，sample1_unpair_clean是生成的去除了一條read的不成對的文件
#ILLUMINACLIP接頭的去除鲸拥，先說明位置拐格，在說明怎么去（2:30:10:1:true）
#SLIDINGWINDOW:4:20 滑窗設(shè)置的是4bp，如果低于20則被去除掉
#MINLEN:50 將read小于50的去除
#TOPHRED33 將reads的堿基質(zhì)量體系轉(zhuǎn)為phred-33
trimmomatic PE -threads 4 \
sample1_R1.fastq.gz sample1_R2.fastq.gz \
../02clean_data/sample1_paired_clean_R1.fastq.gz \
../02clean_data/sample1_unpair_clean_R1.fastq.gz \
../02clean_data/sample1_paired_clean_R2.fastq.gz \
../02clean_data/sample1_unpair_clean_R2.fastq.gz \
ILLUMINACLIP:/root/miniconda3/share/trimmomatic-0.39-1/adapters/TruSeq3-PE-2.fa:2:30:10:1:true \
LEADING:3 TRAILING:3 \
SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33

80%的reads被留下

image.png

4 序列比對

基于基因組和基于轉(zhuǎn)錄本的序列比對
轉(zhuǎn)錄組和基因組比對區(qū)別

image.png

image.png

新建文件夾 mkdir arab_STAR_genome保存索引文件
先解壓00ref里面的基因組和注釋文件用gunzip刑赶，然后建立索引
整個操作是在rnaseq文件夾下完成

建立索引

# 工作模式是genomeGenerate捏浊，arab_STAR_genome生成的索引文件的位置
#genomeFastaFiles和sjdbGTFfile分別表示參考基因組和注釋信息位置
#sjdbOverhang是可變剪切的預(yù)處理，默認(rèn)值是100撞叨。也可以用read長度減1的方式來做金踪，我們的reads是150.
STAR --runThreadN 6 --runMode genomeGenerate \
--genomeDir arab_STAR_genome \
--genomeFastaFiles 00ref/TAIR10_Chr.all.fasta \
--sjdbGTFfile 00ref/Araport11_GFF3_genes_transposons.201606.gtf \
--sjdbOverhang 149

STAR比對

#genomeDir索引的位置，zcat對所有g(shù)z文件解壓谒所，readFilesIn數(shù)據(jù)清洗后要讀入的進(jìn)行比對數(shù)據(jù)
#outFileNamePrefix輸出位置，并給輸出文件加上sample2_ 前綴
#outSAMtype BAM SortedByCoordinate將SAM格式轉(zhuǎn)換為二進(jìn)制的BAM格式沛申，并進(jìn)行排序
#outBAMsortingThreadN 5 排序采用5線程
#quantMode TranscriptomeSAM GeneCounts將基因組比對結(jié)果轉(zhuǎn)換為轉(zhuǎn)錄組結(jié)果劣领，并計數(shù)每個基因的reads數(shù)。
#TranscriptomeSAM 為了使用RSEM進(jìn)行定量分析做準(zhǔn)備铁材。
STAR --runThreadN 5 --genomeDir arab_STAR_genome \
--readFilesCommand zcat \
--readFilesIn 02clean_data/sample2_paired_clean_R1.fastq.gz \
02clean_data/sample2_paired_clean_R2.fastq.gz \
--outFileNamePrefix 03align_out/sample2_ \
--outSAMtype BAM SortedByCoordinate \
--outBAMsortingThreadN 5 \
--quantMode TranscriptomeSAM GeneCounts

輸出比對bam文件

5 表達(dá)定量

image.png

兩種方式都需要先比對再定量

5.1 STAR+RSEM進(jìn)行定量

 mkdir arab_RSEM   #在rnaseq下新建文件夾arab_RSEM保存resm參考文件 

conda install RSEM  #先安裝RSEM,需要時間比較久
#rsem prepare reference
#需要基因組和注釋文件尖淘，提取參考基因組中每一個轉(zhuǎn)錄組的fasta信息
#arab_RSEM/arab_rsem 保存在arab_RSEM奕锌，加arab_rsem前綴。
rsem-prepare-reference --gtf 00ref/Araport11_GFF3_genes_transposons.201606.gtf \
00ref/TAIR10_Chr.all.fasta \
arab_RSEM/arab_rsem

mkdir 04rsem_out#新建RSEM的輸出文件

#表達(dá)定量
#因為數(shù)據(jù)是雙端測序的所以用paired-end村生，不需要輸出bam文件惊暴，所以no-bam-output
#-p 5線程是5，-q 是輸入bam文件的保存位置
#arab_RSEM/arab_rsem \之前準(zhǔn)備工作的索引位置
#輸出到04rsem_out趁桃，加sample2_rsem前綴
rsem-calculate-expression --paired-end --no-bam-output \
--alignments -p 5 \
-q 03align_out/sample2_Aligned.toTranscriptome.out.bam \
arab_RSEM/arab_rsem \
04rsem_out/sample2_rsem

less sample2_rsem.genes.results |head #查看基于基因的定量結(jié)果

image.png

5.2 STAR+featurecounts進(jìn)行定量
可以STAR后辽话，直接用featurecounts進(jìn)行定量，不做RSEM卫病，將得到的out_counts.txt文件下載到本地油啤，然后進(jìn)行后續(xù)處理

image.png

就在STAR后的結(jié)果文件03align_out中運(yùn)行

conda install subread #安裝featurecounts

#-p就是雙端測序，-a后面跟基因組文件蟀苛，-o后面跟輸出的文件
#-T是線程益咬，-t根據(jù)外顯子進(jìn)行定量，-g輸出gene
#sample*_Aligned.sortedByCoord.out.bam 對當(dāng)前文件中所有sample開頭和bam結(jié)尾的文件進(jìn)行定量
featureCounts -p -a ../00ref/Araport11_GFF3_genes_transposons.201606.gtf -o out_counts.txt -T 6 -t exon -g gene_id sample*_Aligned.sortedByCoord.out.bam

out_counts.txt文件中包含了counts數(shù)

out_counts.txt

6 將RSEM定量后的結(jié)果轉(zhuǎn)換為表達(dá)矩陣

image.png

對RSEM定量分析后得到的sample1_rsem.genes.results 文件進(jìn)行合并帜平，構(gòu)建表達(dá)矩陣

image.png

#對所有當(dāng)前目錄下_rsem.genes.results結(jié)尾的文件進(jìn)行合并幽告，輸出到output.matrix
rsem-generate-data-matrix *_rsem.genes.results >output.matrix
head output.matrix #查看
 wc -l output.matrix #查看文件行數(shù)