?? 注：此文基本全部按照 簡書 劉小澤：轉(zhuǎn)錄組那些事兒 Part II 進(jìn)行赘来，感謝??，以下代碼親測(cè)有效汁汗，如有問題歡迎隨時(shí)與我溝通棋电。
準(zhǔn)備工作??
1. 登錄服務(wù)器（本小白用的是2核8G內(nèi)存的云服務(wù)器）或在自己電腦上操作，下載conda（生信分析下載miniconda3就行）椎组，具體參考linux環(huán)境下的軟件安裝
cd biosoft #進(jìn)入目錄
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)對(duì)應(yī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
【注：鏡像配置時(shí)油狂，如果用的是國外服務(wù)器，直接按下面的代碼添加國際鏡像即可寸癌。如果不添加bioconda channel专筷，很多生信分析的軟件下載時(shí)找不到≌粑】??
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)境變量時(shí)才能使用這些軟件
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}
#同時(shí)創(chuàng)建多個(gè)平行及層級(jí)目錄
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}
【為了避免每次登錄服務(wù)器時(shí)都要重新定義變量味咳，可以將以上變量保存在shell腳本中庇勃，登錄時(shí)激活一下即可??】
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(項(xiàng)目)—>SRS(樣本)—>SRX(數(shù)據(jù)產(chǎn)生)—>SRR(數(shù)據(jù)本身)
具體參考 簡書 劉小澤：轉(zhuǎn)錄組那些事兒 Part II
-
本次選擇的示例數(shù)據(jù)：GSE69175
SRR2038215-SRR2038216: 對(duì)照組
SRR2038217-SRR2038218: 實(shí)驗(yàn)組 - 數(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編號(hào)`一個(gè)一個(gè)下載
??但是測(cè)試發(fā)現(xiàn)國內(nèi)服務(wù)器下載速度非常慢责嚷，國外服務(wù)器可以達(dá)到幾十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 設(shè)置時(shí)間掂铐，避免超時(shí)自動(dòng)停止
#下載速度很慢罕拂，幾kb/s，可以試試本地下載或從國外服務(wù)器下載）
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
??測(cè)試發(fā)現(xiàn)會(huì)報(bào)錯(cuò)： 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文件轉(zhuǎn)為fastq.gz文件**
fastq-dump --gzip --split-e *.sra #sra轉(zhuǎn)為fastq.gz
#其中split-e表示如果是單端測(cè)序則一個(gè)sra文件出來一個(gè)fastq文件烤镐，
#如果是雙末端，則一個(gè)sra文件對(duì)應(yīng)兩個(gè)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各自的合集再整合到一起棍鳖，形成一個(gè)數(shù)據(jù)配置文件
cp raw_conf $CLEAN
注：pfastq-dump據(jù)說比fastq-dump更快炮叶，具體方法參考
1. 簡書：如何進(jìn)行SRA到fastq格式的快速轉(zhuǎn)換
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
了解更多??簡書：使用fastp進(jìn)行數(shù)據(jù)質(zhì)控
4. 下載參考基因組和注釋文件并構(gòu)建索引
從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
#構(gòu)建索引
hisat2-build -p 2 $GENOME/hg19.fa hg19
#p代表線程數(shù)，如果服務(wù)器核數(shù)或內(nèi)存較大可增加線程數(shù)
#運(yùn)行時(shí)間較長渡处，約幾個(gè)小時(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. 比對(duì)
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的三件套：轉(zhuǎn)換(view)镜悉、排序(sort)、建索引(index)
轉(zhuǎn)換: -S指輸入文件格式（不加-S默認(rèn)輸入是bam）医瘫，-b指定輸出文件（默認(rèn)輸出sam）【如果要bam轉(zhuǎn)sam侣肄，-h設(shè)置輸出sam時(shí)帶上頭注釋信息】
排序： 對(duì)bam排序，-@ 線程數(shù)醇份， -o輸出文件
索引： 結(jié)果會(huì)產(chǎn)生.bai文件【必須排序后才能建索引～就像體育課必須從高到矮排好以后再報(bào)數(shù)】
6. 基本信息統(tǒng)計(jì)
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號(hào)染色體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)計(jì)
#先用bam_stat.py對(duì)bam文件統(tǒng)計(jì)，看下比對(duì)率
bam_stat.py -i $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.bam.stat
具體運(yùn)行結(jié)果見 簡書 劉小澤：轉(zhuǎn)錄組那些事兒 Part II
7. reads計(jì)數(shù)
基于基因組水平僚纷，可以用Htseq-count和featureCounts
1）Htseq-count：它是用python寫的一個(gè)腳本矩距，conda install -c bcbio htseq -y安裝好以后可以直接拿來用【運(yùn)行約幾十分鐘】
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）對(duì)兩個(gè)軟件的結(jié)果進(jìn)行合并
##合并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)計(jì)一下兩個(gè)軟件的計(jì)數(shù)之和
#統(tǒng)計(jì)featureCounts
perl -alne '$sum += $F[2]; END {print $sum}' matrix_2.count
#結(jié)果是1880017
#統(tǒng)計(jì)htseq-count,結(jié)果是2863201
#我的統(tǒng)計(jì)結(jié)果與原文有些差別怖竭，不知是否由于軟件安裝版本不同導(dǎo)致
具體參數(shù)描述及運(yùn)行結(jié)果見 簡書 劉小澤：轉(zhuǎn)錄組那些事兒 Part II
在我的不懈努力（折騰了我的Mac加兩個(gè)服務(wù)器??）下锥债，目前基本流程能夠運(yùn)行下來，
下一步：
- 詳細(xì)了解數(shù)據(jù)背后的含義痊臭；
- 差異基因的篩選及用R進(jìn)行可視化哮肚。
??～～寫在最后的一些不相關(guān)～～??
1?? 關(guān)于軟件的安裝與卸載：
如果直接運(yùn)行了shell腳本，如conda广匙，一般無需更改環(huán)境變量允趟；如果是一般軟件的安裝，如SRA Toolkit則需要自己添加環(huán)境變量：vim ~/.bash_profile 進(jìn)入編輯環(huán)境變量艇潭，在PATH后面添加冒號(hào)加絕對(duì)路徑（一般加到bin文件）拼窥，如:/Users/chenxiaoxi/miniconda3/bin戏蔑，然后source ~/.bash_profile 激活環(huán)境變量。如果卸載SRA Toolkit則需去掉PATH同時(shí)刪除文件鲁纠。
2?? 常用的一些小命令
echo
du -sh *
du -sh
du -sh 文件名
history
history | grep prefetch
3?? 最近學(xué)到的不同服務(wù)器切換以及數(shù)據(jù)遷移的小命令
su chenxi
exit
scp -r hg19.fa.gz chenxi@521:~/ # 將當(dāng)前服務(wù)器的hg19.fa.gz文件遷移至IP地址為521用戶名為chenxi的家目錄下
4?? 一些感悟：
命令的三重點(diǎn)：input总棵、output、process改含；
多用Tab鍵補(bǔ)齊的方式情龄；
時(shí)刻想著自己在哪里。捍壤。骤视。（當(dāng)前目錄）；
多用： 命令-h 或 命令--help 或 man 命令鹃觉；
不知道某個(gè)命令的含義時(shí)就搜索或試著運(yùn)行看看
...
?? ??最后的最后专酗，感謝我的“人肉搜索引擎”小徐同學(xué)非常耐心（幾乎抓狂）的指導(dǎo)????
