使用sra toolkit,hisat2,stringtie进行有参考基因组的拼接

本文以SRR7663112为例，分为以下几步
建立工作目录 $ mkdir 221 && cd 221，目录下有五个子目录分别为
a)$ mkdir genes ——存放基因注释 *.gff
b) $ mkdir genome ——存放完整基因 *.fna
c)$ mkdir indexes ——存放索引 *.ht2 和用来抽取索引的extract_exons.py，extract_splice_sites.py
d) $ mkdir samples ——存放自己测序的样本*_1.fastq.gz和*_2.fastq.gz
e)$ mkdir sam ——存放回贴后的文件 SRR7663112_chrX.sam

1.下载序列 SRR7663112.sra

$ prefetch SRR7663112

这个命令会将文件下载到~//ncbi/public/sra 目录下
可以使用mv命令将文件移动到工作目录下
mv ~ncbi/public/sra/SRR7663112.sra ./samples

2.将双端测序数据分开,使用hisat2将SRR7663112.sra文件转换为
SRR7663112_1.fastaq.gz和SRR7663112_2.fastaq.gz

$ cd samples
$ ls *.sra | xargs fastq-dump --split-3 --gzip
#这个命令将目录下所有.sra 文件转换

脚本 ~/221/samples/fastq-dump.sh

#!/bin/bash
#PBS -N fast1-dum
#PBS -l nodes=1:ppn=4
#PBS -l walltime=5:00:00

cd /home/201199800009/221/samples

ls *.sra | xargs fastq-dump --split-3 --gzip
########################################################
$ qsub fastq-dump.sh

3.使用hisat2构建 index

$ cd indexes
$ extract_splice_sites.py ../genes/chrX.gtf >chrX.ss 
$ extract_exons.py ../genes/chrX.gtf >chrX.exon
$ hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran

脚本~/221/build_indexes.sh

#!/bin/bash
#PBS -N build_index
#PBS -q cu
#PBS -l nodes=1:ppn=4
#PBS -l walltime=2:00:00

cd /home/201199800009/221/indexes

hisat2-build --ss chrX.ss --exon chrX.exon ../genome/GCA_001185155.1_Zosma_marina.v.2.1_genomic.fna chrX_tran

########################################################
$ qsub build_indexes.sh

3.使用hisat2进行mapping,输出sam 文件

$ cd /home/201199800009/221/sam
$ hisat2 -p 8 --dta -x ../indexes/chrX_tran -1 ../samples/SRR7663112_1.fastq.gz -2../samples/SRR7663112_2.fastq.gz -S SRR7663112_chrX.sam 

脚本/home/201199800009/221/sam/mapping.sh

#!/bin/bash
#PBS -N mapping
#PBS -l nodes=1:ppn=4
#PBS -l walltime=999:00:00

cd /home/201199800009/221/sam

hisat2 -p 8 --dta -x ../indexes/chrX_tran -1 ../samples/SRR7663112_1.fastq.gz -2 ../samples/SRR7663112_2.fastq.gz -S SRR7663112_chrX.sam 
##########################################################################
$ qsub mapping.sh

4.使用samtools把sam文件转换为bam文件

samtools sort -@ 8 -o  SRR7663112_chrX.bam SRR7663112_chrX.sam

脚本/home/201199800009/221/sam/sam_to_bam.sh

#!/bin/bash
#PBS -N sam_to_bam
#PBS -l nodes=1:ppn=4
#PBS -l walltime=999:00:00

cd /home/201199800009/221/sam

samtools sort -@ 8 -o  SRR7663113_chrX.bam SRR7663113_chrX.sam
samtools sort -@ 8 -o  SRR7663114_chrX.bam SRR7663114_chrX.sam
samtools sort -@ 8 -o  SRR7663112_chrX.bam SRR7663112_chrX.sam
################################
qsub sam/sam_to_bam.sh

5.使用stringtie 进行拼接和表达水平评价

$ stringtie -p 8 -G ./genes/GCA_001185155.1_Zosma_marina.v.2.1_genomic.gff -o SRR7663113_chrX.gtf -l SRR7663113 ./sam/SRR7663113_chrX.bam

输出 SRR7663112_chrX.gtf 

6.将所有转录本合并（merge)，merge步骤可以创建出所有样本里面都有的转录本，方便下一步的对比
首先创建一个文件vim mergelist.txt内容是需要合并的所有的文件名

 SRR7663112_chrX.gtf
 SRR7663113_chrX.gtf
 SRR7663114_chrX.gtf

 $ stringtie --merge -p 8 -G ./genes/ GCA_001185155.1_Zosma_marina.v.2.1_genomic.gff -o stringtie_merged.gtf  ./mergelist.txt

7.重新组装转录本并估算基因表达丰度，并为ballgown创建读入文件

$ stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/SRR7663112_chrX/SRR7663112_chrX.gtf ./sam/SSR7663112_chrX.bam
$ stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/SRR7663113_chrX/SRR7663113_chrX.gtf ./sam/SSR7663113_chrX.bam
$ stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/SRR7663114_chrX/SRR7663114_chrX.gtf ./sam/SSR7663114_chrX.bam