06-RNA-seq数据分析全部流程-从拿到原始测序数据开始，从上游到下游分析，看这一篇就够了

1.准备工作

先建好文件夹，这一部分非常非常非常重要，拥有一个优秀的工作习惯比什么都重要

mkdir 293T-6-RNA-seq
cd 293T-6-RNA-seq
mkdir -p data/rawdata data/cleandata/trim_galore data/cleandata/fastp
mkdir -p Expression/featureCounts Expression/Salmon
mkdir Diff_analysis
mkdir -p Mapping/Hisat2 Mapping/Subjunc
tree

2.数据质量评估（拿到原始数据）

1.fastqc

目标：使用fastqc对原始数据进行质量评估
如果想查看历史命令记录，可以使用history | less或者history | tail

# 激活conda环境
conda env list
conda activate RNA_Seq
# 进入到到自己的文件夹，使用FastQC软件对fastq文件进行质量评估，结果输出到qc/文件夹下
vim qc.sh
i
fastqc -t 6 -o ./ *.fq.gz
# 按esc
:wq
nohup bash qc.sh >qc.log &
# 报告整合
multiqc *zip

cat qc.sh
cat qc.log
# 质控的结果给每个fq文件生成一个zip文件和一个html网页报告
# zip文件里是一些报告图片等，主要需要看的是html网页报告，可以将其下载到本地通过浏览器打开

3.数据过滤

如何判断你的数据需要过滤？

原始序列质量控制的标准为：
去除含接头的reads
过滤去除低质量值的数据，确保数据质量
去除含有N（无法确定碱基信息）的比例大于5%的reads

trim_galore过滤

# 激活小环境
conda activate rna
# 到trim_galore文件夹下先生成一个变量,为样本ID
ls /Data/wu_lab/zyy/project/Human-16-Asthma-Trans/data/rawdata/*_1.fq.gz |awk -F '/' '{print $NF}' | cut -d '_' -f 1,2 >ID

cat ID

#多个样本 vim trim_galore.sh，以下为sh的内容
将命令写入sh脚本，使用nohup &运行sh脚本，适用比较长的脚本

vim trim_galore.sh

i
rawdata=/Data/wu_lab/zyy/project/Human-16-Asthma-Trans/data/rawdata
cleandata=/Data/wu_lab/zyy/project/Human-16-Asthma-Trans/data/cleandata/trim_galore
ID=/Data/wu_lab/zyy/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/ID
cat $ID | while read id
do
  trim_galore -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ${cleandata}    ${rawdata}/${id}_1.fq.gz     ${rawdata}/${id}_2.fq.gz
done

# 按esc
:wq

# 提交任务到后台 可以 用bash或者sh都行 
nohup bash trim_galore.sh >trim_galore.log &

# 如何知道运行完了？
cat trim_galore.log

# 可以用jobs查看，也可以用htop查看

# 使用MultiQc整合FastQC结果
multiqc *.zip

fastp过滤

#先到fastp文件夹下
vim fastp.sh
cat ../trim_galore/ID | while read id
do
fastp -l 20 -q 20 --compression=6 \
  -i ${rawdata}/${id}_R1.fq.gz \
  -I ${rawdata}/${id}_R2.fq.gz \
  -o ${cleandata}/${id}_clean_1.fq.gz \
  -O ${cleandata}/${id}_clean_2.fq.gz \
  -R ${cleandata}/${id} \
  -h ${cleandata}/${id}.fastp.html \
  -j ${cleandata}/${id}.fastp.json
done
# 运行fastp脚本
nohup bash fastp.sh >fastp.log &
# fastp同样会生成html文件，下载到本地进行查看

# 整合成一个文件
multiqc *json

数据过滤前后比较

#过滤前后reads数不会变
# 进入过滤目录
cd /Data/wu_lab/zyy/project/293T-shF10-6-RNA-seq/data/cleandata/trim_galore/

# 原始数据
##查看过滤后的reads数
zcat /Data/wu_lab/zyy/project/293T-shF10-6-RNA-seq/data/rawdata/293-0-1_R1.fq.gz |paste - - - - | wc -l
zcat /Data/wu_lab/zyy/project/293T-shF10-6-RNA-seq/data/rawdata/293-0-1_R2.fq.gz |paste - - - - | wc -l
##过滤后的R1和R2的剪辑个数不一样了
zcat 293-0-1_R1_val_1.fq.gz |paste - - - - |cut -f 2 |tr -d '\n' |wc -m

zcat 293-0-1_R2_val_2.fq.gz |paste - - - - |cut -f 2 |tr -d '\n' |wc -m

#  过滤前的数据
zcat ../../rawdata/293-0-1_R1.fq.gz |paste - - - - > raw.txt

#  过滤后的数据
zcat 293-0-1_R1_val_1.fq.gz |paste - - - - > trim.txt
awk '(length($4)<63){print$1}' trim.txt > Seq.ID

head -n 100 Seq.ID > ID100
grep -w -f ID100 trim.txt | awk '{print$1,$4}' > trim.sm
grep -w -f ID100 raw.txt | awk '{print$1,$4}' > raw.sm
paste raw.sm trim.sm | awk '{print$2,$4}' | tr ' ' '\n' |less -S

4.数据比对

目标：使用两个软件对fq数据进行比对，得到比对文件sam/bam，并探索比对结果。

需要准备：
1.参考基因组文件fasta
2.参考基因组注释文件gff/gtf

常用参考基因组数据库
参考基因组准备:注意参考基因组版本信息，可以用ncbi或者Ensembl数据库，一般用Ensembl数据库，更新较快，基因没有重名

Ensembl
添加链接描述

NCBI
添加链接描述

UCSC
添加链接描述

4.1参考基因组

## 参考基因组准备:注意参考基因组版本信息
# 进入到参考基因组目录
cd /Data/wu_lab/zyy/database/GRCh38.p14/

# 下载基因组序列axel  curl  
nohup axel -n 100 https://ftp.ensembl.org/pub/release-113/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz > dna.log &

# 下载转录组序列
nohup axel -n 100  http://ftp.ensembl.org/pub/release-111/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz >rna.log &

# 下载基因组注释文件
nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/gtf/homo_sapiens/Homo_sapiens.GRCh38.111.gtf.gz >gtf.log &

nohup axel -n 100 https://ftp.ensembl.org/pub/release-111/gff3/homo_sapiens/Homo_sapiens.GRCh38.111.gff3.gz >gff.log&

4.2数据比对的过程

1.建索引：为了将短片段快速比对到基因组上的某一个位置
2.比对参考基因组，结果生成sam文件
3.sam转bam
4.bam建索引

4.2.1比对:hisat2

## ----1.构建索引
# 进入参考基因组目录
cd $HOME/database/GRCh38.105
# Hisat2构建索引，构建索引时间比较长，建议提交后台运行，一般会运行20多分钟左右
## 后续索引可直接使用服务器上已经构建好的进行练习
# vim Hisat2Index.sh
mkdir Hisat2Index
fa=Homo_sapiens.GRCh38.dna.primary_assembly.fa
fa_baseName=GRCh38.dna
hisat2-build -p 5 -f ${fa} Hisat2Index/${fa_baseName}

# 运行
nohup bash Hisat2Index.sh >Hisat2Index.sh.log &

## ----2.比对
# 进入比对文件夹
cd /Data/wu_lab/zyy/project/Human-16-Asthma-Trans/Mapping/Hisat2/

## 单个样本比对，步骤分解
index=/Data/wu_lab/zyy/database/GRCm39/Hisat2Index/GRC39.dna
inputdir=/Data/wu_lab/zyy/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/ 
outdir=/Data/wu_lab/zyy/project/Human-16-Asthma-Trans/Mapping/Hisat2 

hisat2 -p 5  -x  ${index} \    
-1 ${inputdir}/AJV_35_1_val_1.fq.gz  \      
-2 ${inputdir}/AJV_35_2_val_2.fq.gz  \      
-S ${outdir}/AJV_35.Hisat_aln.sam


## ----3.sam转bam
samtools sort --threads 5 -o AJV_35.Hisat_aln.sorted.bam AJV_35.Hisat_aln.sam

## ----4.对bam建索引
samtools index AJV_35.Hisat_aln.sorted.bam AJV_35.Hisat_aln.sorted.bam.bai

# 多个样本批量进行比对，排序，建索引
# Hisat.sh内容： 注意命令中的-，表示占位符，表示|管道符前面的输出。
## 此处索引直接使用服务器上已经构建好的进行练习
##1.比对参考基因组；2.sam转bam;3.bam建索引，三个任务串联
# vim Hisat.sh
index=/Data/wu_lab/zyy/database/GRCm39/Hisat2Index/GRC39.dna
inputdir=/Data/wu_lab/zyy/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/ 
outdir=/Data/wu_lab/zyy/project/Human-16-Asthma-Trans/Mapping/Hisat2 

cat ../../data/cleandata/trim_galore/ID | while read id
do
hisat2 -p 5 -x ${index} -1 ${inputdir}/${id}_1_val_1.fq.gz -2 ${inputdir}/${id}_2_val_2.fq.gz 2>${id}.log  | samtools sort -@ 3 -o ${outdir}/${id}.Hisat_aln.sorted.bam - &&  samtools index ${outdir}/${id}.Hisat_aln.sorted.bam ${outdir}/${id}.Hisat_aln.sorted.bam.bai
done
# 统计比对情况
multiqc -o ./ *.log

# 提交后台运行
nohup bash Hisat.sh >Hisat.log &
#查看运行
cat Hisat.log

##也可以