自己在參考大牛的基礎(chǔ)上荸恕，總結(jié)了RNAseq比對流程pipeline腳本。在做好前期準備工作后艰赞，可以一步完成從fastq到表達矩陣的所有步驟莱褒，目前僅支持human Bulk RNAseq數(shù)據(jù)，后續(xù)有機會會繼續(xù)學(xué)習(xí)荒叼、完善轿偎，擴展更多選項。十分歡迎大家參考使用被廓，并提供建議與意見坏晦。
1、準備分析環(huán)境
1.1 linux文件夾環(huán)境
git clone https://github.com/lishensuo/fq2count.git
# gitee備用(快)： git clone https://gitee.com/li-shensuo/fq2count.git
tar -xvf fq2count.tar
work_path=/home/ssli/fq2count
cd $work_path
#為script腳本增加執(zhí)行權(quán)限
##前一個為批量比對的腳本文件伊者；后四個為去除rRNA相關(guān)的腳本文件
chmod u+x ${work_path}/scripts/BulkRNAseq.sh
chmod u+x ${work_path}/scripts/indexdb_rna
chmod u+x ${work_path}/scripts/merge-paired-reads.sh
chmod u+x ${work_path}/scripts/sortmerna
chmod u+x ${work_path}/scripts/unmerge-paired-reads.sh
#新建將會用到的文件夾
mkdir ${work_path}/0.rawfq
mkdir ${work_path}/1.rm_rrna
mkdir ${work_path}/2.trim
mkdir ${work_path}/3.align
mkdir ${work_path}/4.count]
mkdir ${work_path}/rrna_index
1.2 conda環(huán)境
cat requirement.txt
# bioconda::samtools==1.14
# bioconda::subread==2.0.1
# bioconda::refgenie==0.12.1
# bioconda::sra-tools==2.11.0
# bioconda::hisat2==2.2.1
conda create -n fq2count -y
conda activate fq2count
conda install -c conda-forge mamba -y
conda install --file=requirement.txt -y
2英遭、準備數(shù)據(jù)
2.1 參考基因組gtf文件
- 通過refgenie下載，這里以hg38版本為例
#第一次使用refgenie需要運行下面兩行代碼
mkdir ~/refgenie
refgenie init -c ~/refgenie/genome_config.yaml
#下載hg38 gtf
refgenie pull hg38/gencode_gtf -c ~/refgenie/genome_config.yaml
echo $(refgenie seek hg38/gencode_gtf -c ~/refgenie/genome_config.yaml)
2.2 比對軟件(hisat2)索引文件
- 通過refgenie下載亦渗，這里以hg38版本為例
refgenie pull hg38/hisat2_index -c ~/refgenie/genome_config.yaml
echo $(refgenie seek hg38/hisat2_index -c ~/refgenie/genome_config.yaml)
2.3 構(gòu)建rRNA索引文件（optional）
- 如果考慮去除fastq文件里的rRNA reads挖诸，需要運行這一步
#使用腳本文件scripts/indexdb_rna創(chuàng)建核糖體索引
cd ${work_path}/rrna_index
${work_path}/scripts/indexdb_rna --ref ${work_path}/scripts/rRNA_databases/silva-bac-16s-id90.fasta,./silva-bac-16s-db:\
${work_path}/scripts/rRNA_databases/silva-bac-23s-id98.fasta,./silva-bac-23s-db:\
${work_path}/scripts/rRNA_databases/silva-arc-16s-id95.fasta,./silva-arc-16s-db:\
${work_path}/scripts/rRNA_databases/silva-arc-23s-id98.fasta,./silva-arc-23s-db:\
${work_path}/scripts/rRNA_databases/silva-euk-18s-id95.fasta,./silva-euk-18s-db:\
${work_path}/scripts/rRNA_databases/silva-euk-28s-id98.fasta,./silva-euk-28s:\
${work_path}/scripts/rRNA_databases/rfam-5s-database-id98.fasta,./rfam-5s-db:\
${work_path}/scripts/rRNA_databases/rfam-5.8s-database-id98.fasta,./rfam-5.8s-db
2.4 測序數(shù)據(jù)
#如下為示例分析文件
cat SRR_list.txt
# SRR12720999
# SRR12721000
# SRR12721001
SRR_file=SRR_list.txt
cd ${work_path}/0.rawfq
for srr in $(cat $SRR_file)
do
#下載sra文件
prefetch -p -X 35G ${srr} -O .
#拆分為fastq文件
fasterq-dump --split-files ${srr} -p -O ./
rm -rf ${srr}
done
#注意是未壓縮的fastq文件，主要是考慮需要去除rRNA的情況
ls
-rw-r--r-- 1 ssli 7.5G Feb 20 10:34 SRR12720999_1.fastq
-rw-r--r-- 1 ssli 7.5G Feb 20 10:34 SRR12720999_2.fastq
-rw-r--r-- 1 ssli 7.4G Feb 20 10:45 SRR12721000_1.fastq
-rw-r--r-- 1 ssli 7.4G Feb 20 10:45 SRR12721000_2.fastq
-rw-r--r-- 1 ssli 8.0G Feb 20 10:53 SRR12721001_1.fastq
-rw-r--r-- 1 ssli 8.0G Feb 20 10:53 SRR12721001_2.fastq
3法精、批量比對
- 指定相關(guān)參數(shù)
work_path=/home/ssli/rnaseq/fq2count
#交代是否去除核糖體rRNA：0表示不去除多律，1表示去除(比較耗時)
Trim_rRNA=0
#交代測序讀長（根據(jù)實際fastq數(shù)據(jù)）
read_length=100
#交代線程數(shù)
threads=20
#交代參考文件的路徑（根據(jù)需要選擇合適的版本）
hisat_index=$(refgenie seek hg38/hisat2_index -c ~/refgenie/genome_config.yaml)
gtf_path=$(refgenie seek hg38/gencode_gtf -c ~/refgenie/genome_config.yaml)
- 批量比對
${work_path}/scripts/BulkRNAseq.sh $work_path $Trim_rRNA $read_length \
$threads $hisat_index $gtf_path
head ${work_path}/4.count/expression_matrix.txt
# Geneid ../3.align/SRR12720999_nsorted.bam ../3.align/SRR12721000_nsorted.bam ../3.align/SRR12721001_nsorted.bam
# ENSG00000223972.5 0 0 0
# ENSG00000227232.5 24 15 14
# ENSG00000278267.1 12 5 8
# ENSG00000243485.5 0 0 1
# ENSG00000284332.1 0 0 0
# ENSG00000237613.2 0 0 0
# ENSG00000268020.3 0 0 0
# ENSG00000240361.2 0 0 0
# ENSG00000186092.6 0 0 0
