轉錄組分析入門 2 —— 基本流程

?? 注：此文基本全部按照 簡書 劉小澤：轉錄組那些事兒 Part II 進行，感謝??，以下代碼親測有效，如有問題歡迎隨時與我溝通。

準備工作??

1. 登錄服務器（本小白用的是2核8G內(nèi)存的云服務器）或在自己電腦上操作，下載conda（生信分析下載miniconda3就行），具體參考linux環(huán)境下的軟件安裝 
cd biosoft #進入目錄
uname -a #查看系統(tǒng)版本
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh #下載conda
bash Miniconda3-latest-Linux-x86_64.sh #安裝系統(tǒng)對應版本的miniconda
source ~/.bashrc #激活conda
#添加清華鏡像
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

【注：鏡像配置時，如果用的是國外服務器，直接按下面的代碼添加國際鏡像即可。如果不添加bioconda channel，很多生信分析的軟件下載時找不到?！??

conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
2. 配置conda，下載軟件
conda create -n rna-seq python=3 samtools fastp sra-tools hisat2 rseqc subread -y 
#創(chuàng)建rna-seq的環(huán)境變量，并下載samtools等軟件，只有激活rna-seq環(huán)境變量時才能使用這些軟件
conda install -c bcbio htseq -y
conda install numpy pysam -y
3. 配置好工作路徑
RNA=/home/chenxi/project/rna-seq/data
mkdir -p $RNA/{raw_data,clean_data,ref/{genome,gtf,index},align,stats,matrix}
#同時創(chuàng)建多個平行及層級目錄
RAW=$RNA/raw_data
CLEAN=$RNA/clean_data
GENOME=$RNA/ref/genome
GTF=$RNA/ref/gtf
INDEX=$RNA/ref/index
ALIGN=$RNA/aign
MATRIX=$RNA/matrix
STATS=$RNA/stats
mkdir -p $MATRIX/{htseq,featureCounts}

【為了避免每次登錄服務器時都要重新定義變量，可以將以上變量保存在shell腳本中，登錄時激活一下即可??】

vim ~/project/rna-seq/env.sh #編輯
source ~/project/rna-seq/env.sh #激活
echo $CLEAN #檢查是否work

分析流程??????

1. 數(shù)據(jù)的下載（從GEO數(shù)據(jù)庫下載SRA原始數(shù)據(jù)）

SRP(項目)—>SRS(樣本)—>SRX(數(shù)據(jù)產(chǎn)生)—>SRR(數(shù)據(jù)本身)
具體參考 簡書 劉小澤：轉錄組那些事兒 Part II

  • 本次選擇的示例數(shù)據(jù)：GSE69175
    SRR2038215-SRR2038216: 對照組
    SRR2038217-SRR2038218: 實驗組
  • 數(shù)據(jù)下載方法

1）NCBI官方的 SRA Toolkit 中的prefetch命令下載

#前面已經(jīng)安裝sra-tools，可以直接用prefetch，如果沒有則需要先去NCBI官網(wǎng)安裝sra-tools
for i in `seq 5 8`;
do
prefetch SRR203821${i}
done
#也可直接`prefetch SRR編號`一個一個下載

??但是測試發(fā)現(xiàn)國內(nèi)服務器下載速度非常慢，國外服務器可以達到幾十Mb/s??

2）aspera 工具下載

wget -T 8000 https://download.asperasoft.com/download/sw/connect/3.9.8/ibm-aspera-connect-3.9.8.176272-linux-g2.12-64.tar.gz 
#下載aspc軟件，-T 設置時間，避免超時自動停止
#下載速度很慢，幾kb/s，可以試試本地下載或從國外服務器下載）
gunzip ibm-aspera-connect-3.9.8.176272-linux-g2.12-64.tar.gz #解壓ascp
source ~/.bash_profile #激活ascp
#使用ascp下載sra數(shù)據(jù)??
for i in `seq 15 18`;
do 
ascp -v -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
-k 1 -T -l200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR203/SRR20382${i}/SRR20382${i}.sra ./ 
done

??測試發(fā)現(xiàn)會報錯： Failed to open TCP connection for SSH, 目前還未找到原因；
但是用ascp下載EBI的數(shù)據(jù)灰常好用??，下載NCBI的數(shù)據(jù)貌似不太好用。

ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR100/070/SRR10099870 ./

3）wget, curl 命令直接下載

了解更多：
下載NCBI SRA數(shù)據(jù)的最佳方法（來自知乎）
簡書：SRA 數(shù)據(jù)下載自救指南
簡書：安裝Aspera Connect工具下載sra數(shù)據(jù)

2. 數(shù)據(jù)下載完成以后用fastq-dump將sra文件轉為fastq.gz文件**
fastq-dump --gzip --split-e *.sra #sra轉為fastq.gz
#其中split-e表示如果是單端測序則一個sra文件出來一個fastq文件，
#如果是雙末端，則一個sra文件對應兩個fastq文SRRXXXXXX_1.fastq,SRRXXXXXX_2.fastq
find $RAW -name "*.gz" | sort | grep 1.fastq.gz >1
find $RAW -name "*.gz" | sort | grep 2.fastq.gz >2
paste 1 2 > raw_conf 
#將read1和read2各自的合集再整合到一起，形成一個數(shù)據(jù)配置文件 
cp raw_conf $CLEAN
fq.gz 文件

注：pfastq-dump據(jù)說比fastq-dump更快，具體方法參考
1. 簡書：如何進行SRA到fastq格式的快速轉換
2. git pfastq-dump

3. 質(zhì)控過濾
cd $CLEAN
cat raw_conf | while read id;
do 
line=(${id}); 
fq1=${line[0]}; fq2=${line[1]}; 
fastp -i $fq1 -I $fq2 -o out.$(basename $fq1) -O out.$(basename $fq2) -w 2; 
done
結果文件 clean_data

了解更多??簡書：使用fastp進行數(shù)據(jù)質(zhì)控

4. 下載參考基因組和注釋文件并構建索引

從UCSC數(shù)據(jù)庫下載參考基因組文件：https://hgdownload.soe.ucsc.edu/downloads.html
從GENCODE下載注釋信息：https://www.gencodegenes.org/

##下載 hg19 基因組（解壓后大小約3G）
cd $GENOME
for i in $(seq 1 22) X Y M;
do
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do
cat chr${i}.fa >> hg19.fa;
done
rm chr*
#或者可以直接下載官網(wǎng)的成品??
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
#下載注釋文件（解壓后大小約1.3G）
cd GTF
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_33/gencode.v33.annotation.gtf.gz
gunzip *.gz
#構建索引
hisat2-build -p 2 $GENOME/hg19.fa hg19
#p代表線程數(shù)，如果服務器核數(shù)或內(nèi)存較大可增加線程數(shù)
#運行時間較長，約幾個小時
#也可以從hisat2官網(wǎng)直接下載索引文件??
wget https://cloud.biohpc.swmed.edu/index.php/s/hg19/download
mv download hg19.tar.gz #文件重命名
tar -zxvf hg19.tar.gz #解壓
索引文件
5. 比對
for i in `seq 15 18`;
do 
hisat2 -t -p 2 -x $INDEX/hg19 \
-1 $CLEAN/out.SRR20382${i}_1.fastq.gz \ 
-2 $CLEAN/out.SRR20382${i}_2.fastq.gz \
-S SRR20382${i}.sam
samtools view -Sb SRR20382${i}.sam > SRR20382${i}.bam
samtools sort -@ 2 -o SRR20382${i}.sorted.bam SRR20382${i}.bam
samtools index SRR20382${i}.sorted.bam; rm *.sam
done
比對后的結果文件

使用了samtools的三件套：轉換(view)、排序(sort)、建索引(index)
轉換: -S指輸入文件格式（不加-S默認輸入是bam），-b指定輸出文件（默認輸出sam）【如果要bam轉sam，-h設置輸出sam時帶上頭注釋信息】
排序： 對bam排序，-@ 線程數(shù)， -o輸出文件
索引： 結果會產(chǎn)生.bai文件【必須排序后才能建索引～就像體育課必須從高到矮排好以后再報數(shù)】

6. 基本信息統(tǒng)計
cd $STATS
for i in `seq 15 18`;
do
samtools flagstat $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.sorted.flag
done
#如果想根據(jù)flag提取特定區(qū)域，比如想查看1號染色體100-10000的區(qū)域的信息
#samtools view -b -f 0x10 $ALIGN/SRR20382${i}.sorted.bam chr1:100-10000 > ${i}.flag.bam
#samtools flagstat ${i}.flag.bam

#使用RSeqQC統(tǒng)計
#先用bam_stat.py對bam文件統(tǒng)計，看下比對率
bam_stat.py -i $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.bam.stat

具體運行結果見 簡書 劉小澤：轉錄組那些事兒 Part II

7. reads計數(shù)

基于基因組水平，可以用Htseq-count和featureCounts

1）Htseq-count：它是用python寫的一個腳本，conda install -c bcbio htseq -y安裝好以后可以直接拿來用【運行約幾十分鐘】

cd $MATRIX/htseq
for i in `seq 15 18`;
do
htseq-count -s no -r name -f bam $ALIGN/SRR20382${i}.sorted.bam \
$GTF/gencode.v33.annotation.gtf \
>SRR20382${i}.count 2>SRR20382${i}.log
done

2）featureCounts：隸屬于subread套件【相比于htseq更快，約幾分鐘】

cd $MATRIX/featureCounts
begin=$(date +%s)
for i in `seq 15 18`;
do 
featureCounts -T 2 -p -t exon -g gene_id -a $GTF/gencode.v33.annotation.gtf \
-o SRR20382${i}.feature.count $ALIGN/SRR20382${i}.bam; 
done
tim=$(echo "$(date +%s)-$begin" | bc)
printf "time of featureCounts for 4 samples: %.4f seconds" $tim

3）對兩個軟件的結果進行合并

##合并htseq生成的count文件到matrix.count
cd $MATRIX/htseq
perl -lne 'if ($ARGV=~/(.*).count/){print "$1\t$_"}' *.count | grep -v "_" >matrix.count
##合并featureCounts生成的count文件到matrix_2.count
cd $MATRIX/featureCounts
for i in `seq 15 18`;do
cut -f 1,7 SRR20382${i}.feature.count | grep -v "^#" > SRR20382${i}.matrix
sed -i '1d' SRR20382${i}.matrix
done
perl -lne 'if ($ARGV=~/(.*).matrix/){print "$1\t$_"}' *.matrix >matrix_2.count

4）統(tǒng)計一下兩個軟件的計數(shù)之和

#統(tǒng)計featureCounts
perl -alne '$sum += $F[2]; END {print $sum}' matrix_2.count
#結果是1880017
#統(tǒng)計htseq-count,結果是2863201
#我的統(tǒng)計結果與原文有些差別，不知是否由于軟件安裝版本不同導致

具體參數(shù)描述及運行結果見 簡書 劉小澤：轉錄組那些事兒 Part II

在我的不懈努力（折騰了我的Mac加兩個服務器??）下，目前基本流程能夠運行下來，
下一步：

  1. 詳細了解數(shù)據(jù)背后的含義；
  2. 差異基因的篩選及用R進行可視化。

??～～寫在最后的一些不相關～～??

1?? 關于軟件的安裝與卸載：
如果直接運行了shell腳本，如conda，一般無需更改環(huán)境變量；如果是一般軟件的安裝，如SRA Toolkit則需要自己添加環(huán)境變量：vim ~/.bash_profile 進入編輯環(huán)境變量，在PATH后面添加冒號加絕對路徑（一般加到bin文件），如:/Users/chenxiaoxi/miniconda3/bin，然后source ~/.bash_profile 激活環(huán)境變量。如果卸載SRA Toolkit則需去掉PATH同時刪除文件。
2?? 常用的一些小命令
echo
du -sh *
du -sh
du -sh 文件名
history
history | grep prefetch
3?? 最近學到的不同服務器切換以及數(shù)據(jù)遷移的小命令
su chenxi
exit
scp -r hg19.fa.gz chenxi@521:~/ # 將當前服務器的hg19.fa.gz文件遷移至IP地址為521用戶名為chenxi的家目錄下
4?? 一些感悟：
命令的三重點：input、output、process；
多用Tab鍵補齊的方式；
時刻想著自己在哪里。。。（當前目錄）；
多用： 命令-h 或 命令--help 或 man 命令；
不知道某個命令的含義時就搜索或試著運行看看
...

?? ??最后的最后，感謝我的“人肉搜索引擎”小徐同學非常耐心（幾乎抓狂）的指導????

最后編輯于 ：
?著作權歸作者所有,轉載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容