转录组——上游分析

### 3.1 fastqc

目标：使用fastqc对原始数据进行质量评估

# 激活conda环境
conda activate rna

# 连接数据到自己的文件夹（如果前面已经链接过，这里就不用链接了）
cd $HOME/project/Human-16-Asthma-Trans/data/rawdata

# 链接
ln -s /home/t_rna/data/airway/fastq_raw25000/*gz  ./

# 解压缩（补充）
gzip -c SRR1039510_1.fastq.gz > SRR1039510_1.fastq

# 查看软件安装的路径（补充）
which fastqc

### 3.2 **运行：nohup 与 & 提交后台**

+   **任务查看**：使用 jobs、ps fxww、htop命令查看任务状态

+   **任务进度**：使用less 命令查看 qc.log 运行动态报错等，目录中是否有正常文件结果生成

# 使用FastQC软件对单个fastq文件进行质量评估，结果输出到qc/文件夹下
nohup fastqc -t 6 -o ./ SRR*.fastq.gz 1>qc.log 2>&1 &

使用MultiQc整合FastQC结果，得到一个所有样本的qc质控 qc_fastqc.html 文件报告，下载到本地使用浏览器打开查看

multiqc *.zip -o ./ -n qc_fastqc 
-o是指输出到哪里，-n是指输出的报告的名字

# 激活小环境
conda activate rna

# 进入到 cleandata 目录
cd $HOME/project/Human-16-Asthma-Trans/data/cleandata/

# trim
# -j：使用cpu数
# -q：低质量Q值
# --length 20：过滤最小read长度，read被剪切掉接头后长度会发生变化，小于这个长度的会被过滤
# --max_n 3：read中N含量最大值，大于这个数的read会被过滤，N默认为 read长度的5%，这里N=63*0.05~=3
# -o：输出目录
# --fastqc：对过滤后的cleandata 再运行一次质量评估 fastqc 
# --stringency：read与接头重叠的碱基个数
trim_galore -j 3 -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ./ ../rawdata/SRR1039510_1.fastq.gz ../rawdata/SRR1039510_2.fastq.gz

# 激活小环境
conda activate rna

# 进入到 cleandata 目录
cd $HOME/project/Human-16-Asthma-Trans/data/cleandata/

# 先生成一个变量,为样本ID
ls ../rawdata/*_1.fastq.gz | awk -F'/' '{print $NF}' | cut -d'_' -f1 >ID
笔记：这里的NF是列数，要加$才能用，而行数NR不用加。
# 查看ID类容
cat ID

# ID内容如下：
# SRR1039510
# SRR1039511
# SRR1039512

# -j:线程数，这里使用6个cpu线程
cat ID | while read id
do
  echo "trim_galore -j 3 -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ./  ../rawdata/${id}_1.fastq.gz ../rawdata/${id}_2.fastq.gz"
done >trim.sh
笔记：原始未缩减数据在 /home/t_rna/data/airway/fastq_raw/

# 安装
conda activate rna
# conda install parafly -y

# 运行
# -CPU：一次投递几个任务，这里一次投递3个样本
nohup ParaFly -c trim.sh -CPU 3 1>trim.log  2>&1 &

# 笔记：
终止任务：kill 
• kill -9 %num：num为jobs命令返回ID
• kill -9 PID：ps fx或者htop可得到PID任务编号)
① 提交任务：nohup bash trim.sh 1>trim.log 2>&1 &
② 使用 ps fx（后面可以接管道符来grep或less该任务；或者 htop 或者 jobs得到任务ID

终止任务：Ctrl+C
只对直接运行的任务有效，后台任务
无效
① 提交任务：bash trim.sh
② 按 ctrl+c

任务暂停：Ctrl+Z
只对直接运行的任务有效，后台任
务无效
① 提交任务：bash trim.sh
② 按 ctrl+z
③ 把暂停的任务放后台：bg %num；把暂停的任务放后台：fg %num 

# 使用MultiQc整合FastQC结果
multiqc *.zip -n qc_trim
# 笔记：查看报告时，过滤前后主要的变化体现在per base N content，我们的序列里N的数目最多只有3个了。

# 进入过滤目录
cd $HOME/project/Human-16-Asthma-Trans/data/cleandata/

# 原始数据
zcat ../rawdata/SRR1039510_1.fastq.gz | paste - - - - > raw.txt

#  过滤后的数据
zcat SRR1039510_1_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

# 笔记：
Linux进阶命令：paste - - - -执行后，四列变成了一行，但用cat -An 还是会显示制表符^I的。
grep -w -f ID100 trim.txt | awk '{print$1,$4}'执行后，这里的awk默认以空格作为分隔符，而不是制表符，所以这里的$4 不是碱基质量那一列。

## 参考基因组准备:注意参考基因组版本信息
# 下载，Ensembl：http://asia.ensembl.org/index.html
# http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/

# 进入到参考基因组目录
mkdir -p $HOME/database/GRCh38.113
cd $HOME/database/GRCh38.113

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

# 下载转录组序列
nohup wget -c https://ftp.ensembl.org/pub/release-113/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz >rna.log &

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

nohup wget -c https://ftp.ensembl.org/pub/release-113/gff3/homo_sapiens/Homo_sapiens.GRCh38.113.gff3.gz >gff.log &

# 上述文件下载完整后，再解压；否则文件不完整就解压会报错
# 再次强调，一定要在文件下载完后再进行解压！！！
# 笔记：看上述文件的下载进展可以用 tail dna.log 这个文件
nohup gunzip *gz >unzip.log &

zless -NS ./project/Human-16-Asthma-Trans/data/rawdata/SRR1039510_1.fastq.gz | paste - - - - | awk '{print $1,$4}' | sed 's/ /\n/g' | sed 's/@/>/g' | less -NS

# 答案1
zless -S SRR1039511_1_val_1.fq.gz |awk '{ if(NR%4==1){print">" substr($0,2)} if(NR%4==2){print}  }' | less -S

# 答案2
zless -S SRR1039510_1_val_1.fq.gz |paste - - - - |cut -f 1,2 |tr '@' '>' |tr '\t' '\n' |less -S

## ----1.构建索引（使用服务器上已经构建好的）
# 进入参考基因组目录
cd $HOME/database/GRCh38.113
# Hisat2构建索引，构建索引时间比较长，建议提交后台运行，一般会运行20多分钟左右
## 后续索引可直接使用服务器上已经构建好的进行练习
# 创建索引存在文件夹
mkdir Hisat2Index
# 笔记：索引占用的空间比较大，大家共用一份就可以了，在/home/t_rna/database/GRCh38.104，参考的基因组文件也在这里面。

# 运行提交后台
nohup hisat2-build -p 5 -f Homo_sapiens.GRCh38.dna.primary_assembly.fa Hisat2Index/GRCh38.dna  1>index.log  2>&1 &

## ----比对
# 进入比对文件夹
cd $HOME/project/Human-16-Asthma-Trans/Mapping/

## 单个样本比对，步骤分解
index=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dna
用echo查看：echo $index
用ls查看：ls $index*

# step1：比对
# -S：输出结果文件，格式为sam格式
# \：表示手动换行，\后面不能有空格
hisat2 -p 5  -x  ${index} \
       -1 ../data/cleandata/SRR1039510_1_val_1.fq.gz \
       -2 ../data/cleandata/SRR1039510_2_val_2.fq.gz \
       -S ./SRR1039510.Hisat_aln.sam

## step2: sam转bam
# -@：使用线程数
samtools sort -@ 5 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam

## step3：对bam建索引，以便后续使用IGV查看
samtools index SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sorted.bam.bai

## ----2.比对
# 进入比对文件夹
cd $HOME/project/Human-16-Asthma-Trans/Mapping/

## 上课使用已经提前构建好的索引，版本为GRCh38.104
## 定义索引前缀index
index=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dna


## 使用echo生成 hisat.sh
cat ../data/cleandata/ID | while read id
do
  echo "hisat2 -p 5 -x ${index} -1 ../data/cleandata/${id}_1_val_1.fq.gz -2 ../data/cleandata/${id}_2_val_2.fq.gz 2>${id}.log |samtools sort -@ 3 -o ${id}.Hisat_aln.sorted.bam - "
done >hisat.sh
笔记：这个sh脚本里面有一个占位符 - ，需要理解一下，也就是用hisat2软件比对后的结果没有输出到sam文件里的这个步骤了（-S ./SRR1039510.Hisat_aln.sam），而是直接传到samtools这个软件来进行分析。
难点：- 占位符 表示前一个任务的输出结果通过管道符传递给后一个命令，并指定位置。
笔记：会随着服务器断开，需要重新赋值定义index，所以生成hisat.sh脚本后，要用cat检查一下命令对不对。

# 运行
# -CPU：一次投递几个任务，这里一次投递3个样本
nohup ParaFly -c hisat.sh -CPU 3 1>hisat.log  2>&1 &

## ----4.对bam建索引
## 使用echo生成 index.sh
cat ../data/cleandata/ID | while read id
do
  echo "samtools index ${id}.Hisat_aln.sorted.bam ${id}.Hisat_aln.sorted.bam.bai"
done >index.sh

# 提交任务到后台 可以 用bash或者sh都行 
nohup ParaFly -c index.sh -CPU 3 1>index.log  2>&1 &

# 统计比对情况
multiqc -o ./ SRR*log -n hisat2.maprate.html

samtools view 查看sam，bam文件；# sam，bam文件转换

cd $HOME/project/Human-16-Asthma-Trans/Expression/

## 定义输入输出文件夹
gtf=/home/t_rna/database/GRCh38.104/Homo_sapiens.GRCh38.104.chr.gtf.gz

# 看一下自己的featureCounts版本，2.0.3 以前的和以后的参数有点改变
featureCounts -v

# featureCounts对bam文件进行计数 v2.0.3以上版本
# -T：线程数
# -p --countReadPairs：这两个参数制定安找双端模式进行定量
# -t exon -g gene_id：按照基因水平进行定量
# -a $gtf：注释文件
# -o all.id.txt：输出的定量表达矩阵

# gene id
nohup featureCounts -T 6 -p --countReadPairs -t exon -g gene_id -a $gtf -o all.id.txt ../Mapping/*.sorted.bam 1>count.log 2>&1 &

# gene name
nohup featureCounts -T 6 -p --countReadPairs -t exon -g gene_name -a $gtf -o all.id_gene_name.txt ../Mapping/*.sorted.bam 1>count_gene_name.log 2>&1 &（代码报错，暂未有解决办法）

# 对定量结果质控
multiqc ./all.id.txt.summary -n qc_count

# 得到表达矩阵
# 处理表头，
less -S all.id.txt |grep -v '#' |cut -f 1,7- |sed 's#../Mapping/##g' |sed 's#.Hisat_aln.sorted.bam##g' >raw_counts.txt

# 使用less -S 
less -S raw_counts.txt

# 也可以 列对齐显示 看看前6行
head raw_counts.txt  |column -t

# 看一下有多少行，去掉表头有多少个基因
wc -l raw_counts.txt

# 看一下有多少列，去掉第一列有多少个样本, awk中NF表示列数，$NF表示最后一列元素
head -n 3 raw_counts.txt| awk -F'\t' '{print NF}'