转录组数据分析-比对

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

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

# 下载基因组序列axel  curl  ##在后台挂上
nohup wget -c http://ftp.ensembl.org/pub/release-105/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz >dna.log &
(rna) Mar402 16:48:59 ~/database/GRCh38.105
$ ll       #查看 dna.log
total 75568
drwxrwxr-x 2 Mar402 Mar402     4096 Apr 23 16:46 ./
drwxrwxr-x 3 Mar402 Mar402     4096 Apr  7 11:12 ../
-rw-rw-r-- 1 Mar402 Mar402     4793 Apr 23 16:50 dna.log
-rw-rw-r-- 1 Mar402 Mar402 77251337 Apr 23 16:50 Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
-rw-rw-r-- 1 Mar402 Mar402   106723 Apr 23 16:50 wget-log
(rna) Mar402 16:50:29 ~/database/GRCh38.105
$ less dna.log  #查看下载进度
(rna) Mar402 16:51:32 ~/database/GRCh38.105
$ ll        #查看 dna.log 大小已经变了
total 99352
drwxrwxr-x 2 Mar402 Mar402      4096 Apr 23 16:46 ./
drwxrwxr-x 3 Mar402 Mar402      4096 Apr  7 11:12 ../
-rw-rw-r-- 1 Mar402 Mar402      7618 Apr 23 16:51 dna.log
-rw-rw-r-- 1 Mar402 Mar402 101580297 Apr 23 16:51 Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
-rw-rw-r-- 1 Mar402 Mar402    139091 Apr 23 16:51 wget-log

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

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

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

# 上述文件下载完整后，再解压；否则文件不完整就解压会报错
# 再次强调，一定要在文件下载完后再进行解压！！！
nohup gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz Homo_sapiens.GRCh38.cdna.all.fa.gz >unzip.log &

# 进入参考基因组目录
cd $HOME/database/GRCh38.105 #进入文件，保证有解压过的fa文件
# 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}   #-p 线程 -f 接的参考序列 接的索引

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

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

## 单个样本比对，步骤分解
index=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dna #输入
inputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/
outdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2 #输出

hisat2 -p 5  -x  ${index} \  #-p 5 5线程 ；-x前缀 ；\手动换行
	   -1 ${inputdir}/SRR1039510_1_val_1.fq.gz \
       -2 ${inputdir}/SRR1039510_2_val_2.fq.gz \
       -S ${outdir}/SRR1039510.Hisat_aln.sam  #匹配的项目文件

samtools sort -@ 5 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam #sort 排序； -@ 线程数；-o生成的bam文件

# vim Hisat.sh
index=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dna
inputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/
outdir=$HOME/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 - 
done#图c

# 统计比对情况
multiqc -o ./ SRR*log #使用multiqc 生成统计图（图E）

# 提交后台运行
nohup sh Hisat.sh >Hisat.log & 
#结果图D

(rna) Mar402 21:10:13 ~/project/Human-16-Asthma-Trans/Mapping/Hisat2
$ samtools view -h SRR1039510.Hisat_aln.sorted.bam | less -S #-h 显示表头

## ----构建索引
# 进入参考基因组目录
cd $HOME/database/GRCh38.105
# subjunc构建索引，构建索引时间比较长大约40分钟左右，建议携程sh脚本提交后台运行
## 后续索引可直接使用服务器上已经构建好的进行练习
# vim SubjuncIndex.sh
mkdir SubjuncIndex
fa=Homo_sapiens.GRCh38.dna.primary_assembly.fa
fa_baseName=GRCh38.dna
subread-buildindex -o SubjuncIndex/${fa_baseName} ${fa}

# 运行
nohup sh SubjuncIndex.sh >SubjuncIndex.sh.log &


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

# vim subjunc.sh
index=/home/t_rna/database/GRCh38.104/SubjuncIndex/GRCh38.dna
inputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore
outdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Subjunc

cat ../../data/cleandata/trim_galore/ID | while read id
do
  subjunc -T 5 -i ${index} -r ${inputdir}/${id}_1_val_1.fq.gz -R ${inputdir}/${id}_2_val_2.fq.gz -o ${outdir}/${id}.Subjunc.bam 1>${outdir}/${id}.Subjunc.log 2>&1 && samtools sort -@ 6 -o ${outdir}/${id}.Subjunc.sorted.bam  ${outdir}/${id}.Subjunc.bam 
done
## &&说明前面的任务和后面的任务是独立的，只有前面的任务运行成功后才运行后面的任务。
# 运行
nohup sh subjunc.sh >subjunc.log &

# 单个样本
samtools flagstat -@ 3 SRR1039510.Hisat_aln.sorted.bam

# 多个样本，vim flagstat.sh
ls *.sorted.bam | while read id
do
  samtools flagstat -@ 3 ${id} > ${id/bam/flagstat}
done
# 质控
multiqc -o ./  *.flagstat

# 运行
nohup sh flagstat.sh >flagstat.log &

##----view查看bam文件
samtools view SRR1039510.Hisat_aln.sorted.bam
samtools view -H SRR1039510.Hisat_aln.sorted.bam
samtools view -h SRR1039510.Hisat_aln.sorted.bam

##----index对bam文件建索引
samtools index SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sorted.bam.bai

##----flagstat统计比对结果
samtools flagstat -@ 3 SRR1039510.Hisat_aln.sorted.bam

##----sort排序 sam转bam并排序
samtools sort -@ 3 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam

##----depth统计测序深度
# 得到的结果中，一共有3列以指标分隔符分隔的数据，第一列为染色体名称，第二列为位点，第三列为覆盖深度
samtools depth SRR1039510.Hisat_aln.sorted.bam >SRR1039510.Hisat_aln.sorted.bam.depth.txt

##----计算某一个基因的测序深度
# 找到基因的坐标
zless -S Homo_sapiens.GRCh38.95.gff3.gz |awk '{if($3=="gene")print}' |grep 'ID=gene:ENSG00000186092' |awk '{print $1"\t"$4"\t"$5}' >ENSG00000186092.bed
samtools depth  -b  ENSG00000186092.bed SRR1039510.Hisat_aln.sorted.bam >ENSG00000186092.bed.depth

# 如何找到多比对的reads，flag值的理解
# (0x100) 代表着多比对情况，所以直接用samtools view -f 0x100可以提取 multiple比对的 情况