2021-07-14:sunruiqi_RNA_seq

今天在新的单位用新的服务器开始处理RNA_seq数据，希望一切顺利啦。

1.建立相应目录

对新数据建立对应实验人员（sunruiqi）、测序类型（RNA_seq）和日期（2021_07_14）的目录。

# 建立后如下：
(base) zexing@bio:~/projects/sunruiqi/RNA_seq/2021_07_14$

# 新建对应的目录
mkdir raw_data clean_data ballgown bam bam_sort sam fastqc_report GSEA MD5_txt  scripts_log

2.使用sratoolkit下载NCBI原始数据并转为FASTA格式

vim新建RNA_seq_script_1将下载原始数据和转换格式放在一起

#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
#     This program is used for RNA-seq data analysis.
# History
#     2021/07/14       zexing            First release

dir=/Data/lizexing/projects/sunruiqi/RNA_seq/2021_07_14/

# 利用for循环进行后续操作
for i in SRR8080937 SRR8080938 SRR8080941 SRR8080942 SRR9302356 SRR9302355 SRR9302354 SRR9302353 SRR9302352 SRR9302351
do

# 对数据使用prefetch下载
prefetch -O ${dir}/clean_data/ ${i}

# 对数据使用fasterq -dump转换格式
此命令未成功，暂时未解决。
fasterq-dump -O ${dir}/clean_data/${i} ${i}/${i}.sra

done

后台运行RNA_seq_script_1：

nohup bash RNA_seq_script_1 > RNA_seq_script_1_log &

3.使用Aspera下载NCBI中相关基因组数据

NCBI查询网址：https://www.ncbi.nlm.nih.gov/genome

使用ascp命令下载相应基因组数据：

# 下载全基因组数据
ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -l 100M -k 1 -T anonftp@ftp.ncbi.nlm.nih.gov:genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/GCF_000003025.6_Sscrofa11.1_genomic.fna.gz ~/reference

使用wget命令下载相应的基因组注释文件：

# 下载基因组注释文件
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/GCF_000003025.6_Sscrofa11.1_genomic.gff.gz

使用hisat2软件建立相应的index

hisat2-build Sscrofa11.1_genomic.fna hisat2_index_Sscrofa11.1

4.在Linux服务器中对RNA_seq数据进行处理

vim新建RNA_seq_script_2将TFAP2C相关数据质控、比对、格式转换、排序、拼接和定量综合在一起。

#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
#     This program is used for RNA-seq data analysis.
# History
#     2021/07/14       zexing            First release
# 设置变量${dir}为常用目录
dir=/Data/lizexing/projects/sunruiqi/RNA_seq/2021_07_14/

# 对数据进行质控
# fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fastq

# 利用for循环进行后续操作
for i in SRR8080937 SRR8080938 SRR8080941 SRR8080942
do
# 对数据进行比对
hisat2 -t -p 16 -x /Data/lizexing/reference/UCSC_hg19/hisat2_index/hisat2_index_hg19 \
-1 ${dir}/clean_data/${i}_1.fastq \
-2 ${dir}/clean_data/${i}_2.fastq \
-S ${dir}/sam/${i}.sam

# 对数据进行格式转换
samtools view -@ 16 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam

# 对数据进行排序
samtools sort -@ 16 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam

# 对数据进行拼接、定量
mkdir ${dir}/ballgown/"$i"
stringtie ${dir}/bam_sort/"$i".bam.sort -o ${dir}/ballgown/"$i"/"$i".gtf \
-p 16 -G /Data/lizexing/reference/UCSC_hg19/hg19_genes.gtf -e -B \
-A ${dir}/ballgown/"$i"/"$i".gene.tab
done

后台运行RNA_seq_script_2：

nohup bash RNA_seq_script_2 > RNA_seq_script_2_log &

vim新建RNA_seq_script_3将PRRSV相关数据质控、比对、格式转换、排序、拼接和定量综合在一起。

#!/bin/bash
# 上面一行宣告这个script的语法使用bash语法，当程序被执行时，能够载入bash的相关环境配置文件。
# Program
#     This program is used for RNA-seq data analysis.
# History
#     2021/07/14       zexing            First release
# 设置变量${dir}为常用目录
dir=/Data/lizexing/projects/sunruiqi/RNA_seq/2021_07_14/

# 对数据进行质控
# fastqc -t 16 -o ${dir}/fastqc_report/ ${dir}/clean_data/*.fastq

# 利用for循环进行后续操作
for i in SRR9302356 SRR9302355 SRR9302354 SRR9302353 SRR9302352 SRR9302351
do
# 对数据进行比对
hisat2 -t -p 16 -x /Data/lizexing/reference/NCBI_sscr_11/hisat2_index/hisat2_index_Sscrofa11.1 \
-1 ${dir}/clean_data/${i}_1.fastq \
-2 ${dir}/clean_data/${i}_2.fastq \
-S ${dir}/sam/${i}.sam

# 对数据进行格式转换
samtools view -@ 16 -S ${dir}/sam/${i}.sam -1b -o ${dir}/bam/${i}.bam

# 对数据进行排序
samtools sort -@ 16 -l 5 -o ${dir}/bam_sort/${i}.bam.sort ${dir}/bam/${i}.bam

# 对数据进行拼接、定量
mkdir ${dir}/ballgown/"$i"
stringtie ${dir}/bam_sort/"$i".bam.sort -o ${dir}/ballgown/"$i"/"$i".gtf \
-p 16 -G /Data/lizexing/reference/NCBI_sscr_11/Sscrofa11.1_genomic.gff -e -B \
-A ${dir}/ballgown/"$i"/"$i".gene.tab

done

后台运行RNA_seq_script_3：

nohup bash RNA_seq_script_3 > RNA_seq_script_3_log &

5.使用prepDE.py脚本提取read_counts数值

• 进入ballgown文件夹，将prepDE.py脚本拷贝至当前文件夹

cp /Data/lizexing/software/prepDE.py3 ./

• 使用python命令直接运行脚本

python prepDE.py3

运行结果中"gene_count_matrix.csv"即是DESeq2的输入文件。