ENA数据库网址：https://www.ebi.ac.uk/ena/browser/view/SRR949627
参考流程1：https://blog.csdn.net/qq_41933915/article/details/81873596
参考流程2：http://www.360doc.com/content/18/0715/20/19913717_770622175.shtml
参考流程3：https://www.jianshu.com/p/ac367d2e275e
参考流程4：https://mp.weixin.qq.com/s?
src=11×tamp=1572415144&ver=1943&signature=oTwCHlt8tEm1vHh5CWPAIRheUw8O3Ejd5zk5IhSaLrSYfZkJgSMp-mIFqk9sAGUAFgDkvIDVGsTMue9xw6XMHE20wJcHErUkg*v1viWn3O3I-DjkTWP6otNJ0lWBryWF&new=1
参考流程5：https://www.jianshu.com/p/970012d1b051
my_env=/mnt/d/WHQ/biotree_example
1、获取sra数据
cd $my_env
mkdir sra
cd sra
prefetch --option-file SRR.list
#注意，此时sra文件保存在prefetch中，可以通过更改配置来更改下载路径，或者等待下载完后复制到sra文件夹内
2、fastq-dump命令获得fastq文件
cd $my_env
mkdir fastq
#将需要处理的文件批量写入一个shell脚本，注意：如果想获得gz文件的话添加参数--gzip
cd fastq
vim fastqdump
#!/bin/bash
for i in `ls ../sra/*.sra`
do
echo "fastq-dump --split-3 -O ./ ${i}"
done
:wq
#将脚本产生的结果重定向至新的shell脚本
./fastqdump > fastqdump_command
#将该结果再后台运行
nohup bash ./fastqdump_command &
#查看是否在运行
jobs
ps
3、使用fastqc进行质检
#安装软件以及java环境
#Java：https://www.oracle.com/technetwork/java/javase/downloads/jre8-downloads-2133155.html
#fastqc：http://www.bioinformatics.babraham.ac.uk/projects/download.html
#multiqc：conda install -c bioconda -c conda-forge multiqc
#解压、配置环境变量、source
cd $my_env
mkdir fastqc
cd ./fastqc
mkdir multiqc
#使用fastqc进行质检 -t 线程数
fastqc -t 10 -o ./ ../fastq/*.fastq#或者是*.fastq.gz
#使用multiqc对fastqc结果进行整合
cd multiqc
multiqc ../*.zip -o ./
#multiqc 使用说明及结果解读：https://www.jianshu.com/p/85da4dcc6020
4、使用trim-galore去接头并质检
cd $my_env
#安装trim-galore
#conda install trim-galore
#使用说明：https://www.jianshu.com/p/7a3de6b8e503
mkdir clean_data
cd clean_data
trim_galore -q 20 --phred33 --stringency 3 --length 20 --gzip -o . `ls ../fastq/SRR*fastq`
4.1、重新做质检和多重质检
cd $my_env
mkdir clean_data_fastqc
cd clean_data_fastqc
fastqc -t 10 -o ./ ../clean_data/*trim*gz
mkdir multiqc
cd multiqc
multiqc ../*zip -o ./
4.2、使用fastp进行质检（作为trim的备选方案）
因为使用trim-galore会导致对read的过于严格的筛选，可能会使hisat2报如下错误：
因此这里使用fastp作为备选方案
微信图片_20191218011729.png
#安装fastp，conda install fastp
cd $my_env
mkdir fastp_trim
cd fastp-trim
fastp -i ../fastq/SRR5125857_1.fastq.gz -o SRR5125857_clean_1.fq.gz -I ../fastq/SRR5125857_2.fastq.gz -O SRR5125857_clean_2.fq.gz
#注意：小写是read1，大写是read2
5、构建索引
cd $my_env
mkdir index
cd index
#下载hisat2 ：conda install hisat2
#获得索引文件
#方法一、使用官网的索引
#人类hg38索引位置：ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/
mkdir prebuild_index-file
cd prebuild_indel_file/
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg38.tar.gz
#方法二、自己下载fa文件和gff3文件，使用hisat2-build建立索引
#总链接：https://www.cnblogs.com/jessepeng/p/9681749.html
#ncbi-genome
#参考基因组：ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_genomic.fna.gz
#注释文件：ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.39_GRCh38.p13/GCF_000001405.39_GRCh38.p13_genomic.gff.gz
#ncbi-assembly
#参考基因组：https://www.ncbi.nlm.nih.gov/projects/r_gencoll/ftp_service/nph-gc-ftp-service.cgi/?HistoryId=NCID_1_40962930_130.14.18.97_5555_1572225083_2630637030_0MetA0_S_HStore&QueryKey=15&ReleaseType=RefSeq&FileType=GENOME_FASTA&Flat=true
#注释文件：https://www.ncbi.nlm.nih.gov/projects/r_gencoll/ftp_service/nph-gc-ftp-service.cgi/?HistoryId=NCID_1_40962930_130.14.18.97_5555_1572225083_2630637030_0MetA0_S_HStore&QueryKey=15&ReleaseType=RefSeq&FileType=GENOME_GFF&Flat=true
#ensembl:帮助链接(https://jingyan.baidu.com/article/4d58d5416548f79dd4e9c036.html)
#总链接：https://asia.ensembl.org/info/data/ftp/index.html
#参考基因组：ftp://ftp.ensembl.org/pub/release-98/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz
#注释文件：ftp://ftp.ensembl.org/pub/release-98/gff3/homo_sapiens/Homo_sapiens.GRCh38.98.gff3.gz
#以enzembl数据库为例
mkdir index_file
cd index_file
wget ftp://ftp.ensembl.org/pub/release-98/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.gz
wget ftp://ftp.ensembl.org/pub/release-98/gff3/homo_sapiens/Homo_sapiens.GRCh38.98.gff3.gz
gunzip Homo_sapiens.GRCh38.98.gff3.gz
gunzip Homo_sapiens.GRCh38.dna.toplevel.fa.gz
#重命名
mv Homo_sapiens.GRCh38.dna.toplevel.fa.gz ref.fa
mv Homo_sapiens.GRCh38.98.gff3.gz ref.gff3
#建立索引
extract_exons_py ref.gff3 >ref.exon
extract_splice_sites.py ref.gff3 >ref.ss
hisat2-build -p 8 ref.fa --ss ref.ss --exon ref.exon ref.fa
#如果内存不够用，可以直接使用下面代码
hisat2-build -p 8 ref.fa ref.fa
6、回帖（mapping）
cd $my_env
mkdir mapping
cd mapping
#单样本处理
hisat2 -p 15 -x ../index/index_file/ref.fa -U ../clean_data/SRR6790713_trimmed.fq.gz -S sample.sam --dta --un-conc-gz sample.unmapped.gz --novel-splicesite-outfile sample.bed >sample.stat
#单端测序批处理
vim hisat2
#!/bin/bash
ls ../clean_data/*trim*gz|while read id
do
echo "hisat2 -p 15 --dta -x ../index/index_file/ref.fa -U ${id} -S ${id##*/}.sam --un-conc-gz ${id##*/}.unmapped.gz --novel-splicesite-outfile ${id##*/}.bed >sample.stat"
done
：wq
./hisat2 >hisat_command
./hisat_command
#注意：单端测序使用-U，双端测序使用-1 -2
#双端测序批处理
#${i%%_*.*} 用于截取获得前缀
vim hisat2
#!/bin/bash
ls ../clean_data/*_1_trimmed.fq.gz|while read id
do
a=${id%%_[0-9]_trimmed.fq.gz}
echo "hisat2 -p 15 --dta -x ../index/index_file/ref.fa -1 ${a}_1_trimmed.fq.gz -2 ${a}_2_trimmed.fq.gz -S ${a##*/}.sam --un-conc-gz ${a##*/}.unmapped.gz --novel-splicesite-outfile ${a##*/}.bed >sample.stat"
done
：wq
./hisat2 >hisat_command
./hisat_command
#${id##*/}中的##*/表示从右边第一个/开始，只保留右边的内容，删掉左边的内容
#注意：--dta是下面使用stringtie的前提，一定要加上！！！！！别落下！！！！
#安装samtools：conda install samtools
mkdir bam
cd bam
#单样本处理
samtools sort -@ 10 --output-fmt BAM -o sample.sort.bam sample.sam
samtools index -@ 5 sample.sort.bam
#批处理
ls ../*.sam | while read id; do echo "samtools sort --output-fmt BAM -o ${id##*/}.sort.bam ${id}"; echo "samtools index ${id##*/}.sort.bam"; done >samtoolscom
./samtoolscom
#注意 ：-@ 线程数
7、拼接和定量stringtie
#安装stringtie：conda install stringtie
#使用说明：https://www.cnblogs.com/wangprince2017/p/9937343.html
#官方应用指南：http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual
cd $my_env
mkdir stringtie
cd stringtie
#少量样本处理
#拼接
stringtie -p 10 -G ../index/index_file/ref.gff3 -o R1.gff3 -l R1 ../mapping/bam/R1.sort.bam
stringtie -p 10 -G ../index/index_file/ref.gff3 -o R2.gff3 -l R2 ../mapping/bam/R2.sort.bam
stringtie -p 10 -G ../index/index_file/ref.gff3 -o R3.gff3 -l R3 ../mapping/bam/R3.sort.bam
stringtie --merge -p 10 -G ../index/index_file/ref.gff3 -o total.merge.gff3 -l merged SRR*.gff3
#定量
stringtie -e -G total.merged.gff3 -o R1.final.gff3 -l R1 ../mapping/R1.sort.bam#不能自己批量定量，只能写脚本定量，参考下面的批量处理
#获得样本列表用于提取矩阵
ls *.final.gff3|perl -ne '@a=split(/\./,$_);print$a[0],"\t",$_;' >gff3.list
#stringtie批量操作
#拼接的批量操作
vim stringtiecom1
#!/bin/bash
ls ../mapping/bam/*.sort.bam|while read id
do
a=${id##*/}
echo "stringtie -p 10 -G ../index/index_file/ref.gff3 -o ${a:0:10}.gff3 -l ${a:0:10} ${id}"
done
./stringtiecom1 >stringtiecom2
./stringtiecom2
#merge操作
stringtie --merge -p 10 -G ../index/index_file/ref.gff3 -o total.merged.gff3 -l merged SRR*.gff3
#定量批量操作
ls ../mapping/bam/*.sort.bam|while read id
do
a=${id##*/}
echo "stringtie -e -G total.merged.gff3 -o ${a:0:10}.final.gff3 -l ${a:0:10} ${id}"
done >stringtiecom3
./stringtiecom3
#获得样本列表用于提取矩阵
ls *.final.gff3|perl -ne '@a=split(/\./,$_);print$a[0],"\t",$_;' >gff3.list
#注意！！！echo中的另一半引号不要落了！
8、导出为基因表达量矩阵和转录本矩阵(未成功)
cd $my_env
mkdir ballgown
cd ballgown
#将7中的结果使用以ballgown识别的结果展示
ls ../mapping/bam/*.sort.bam|while read id
do
a=${id##*/}
echo "mkdir ${a/.sort.bam/} ;cd ${a/.sort.bam/}"
echo "stringtie -B -e -G ../../stringtie/total.merged.gff3 -o ${a:0:10}.final.gtf -l ${a:0:10} ../${id}"
echo "cd .."
done >ballgowncom1
./ballgowncom1
#将数据导入ballgown
setwd(my_env/ballgown)
library(ballgown)
bg<-ballgown(dataDir = "./",samplePattern = "SRR")
#提取表达矩阵,其中有rpkm信息
expr_mat<-texpr(bg,meas="all")
#差异分析（以后做！）
#尝试用python程序提取表达矩阵（未成功，在未merge的时候是可以跑通的）
mkdir ballgown_gtf
cd ballgown_gtf
mv ../*final.gtf .
#使用官网提供的python代码提取基因count矩阵和转录本矩阵
#python代码地址：http://ccb.jhu.edu/software/stringtie/dl/prepDE.py
#将7中创建的sample列表移动到当前路径下,同时把下载的python代码一并移动过来
wget http://ccb.jhu.edu/software/stringtie/dl/prepDE.py
ls *.final.gtf|perl -ne '@a=split(/\./,$_);print$a[0],"\t",$_;' >gtf.list
python2 ./prepDE.py -i ./gtf.list
#报错信息：
#0 SRR8980083
#1 SRR8980084
#Traceback (most recent call last):
# File "./prepDE.py", line 255, in <module>
# geneDict.setdefault(geneIDs[i],{}) #gene_id
#KeyError: 'transcript:ENST00000511562'
8.1 使用featurecount定量获得表达矩阵
mkdir featurecount
cd featurecount
#安装subread：conda install subread 或者是 apt-get install subread
#featurecounts参数解释：http://www.bioinfo-scrounger.com/archives/407
cp ../mapping/bam/*.sort.bam .
featureCounts -T 10 -p -g gene -t exon -a ../index/index_file/ref.gff3 -o ./count *.sort.bam
#提取表达量（也可以不做）
cut -f 1,7 count |grep -v '^#' >feacturCounts.txt#提取第一列（基因名）和第七列（表达量）
#注意：
#1、-g的输入必须在gff文件中有这个内容才可以，默认是gene_id ，但是我下载的没有这个名字，所以用它存在的Name代替
#2、-a表示之前下载的gff文件，stringtie合并后的也能用，但是不知道是否正确
#3、-o表示将定量结果输出到count文件中
#4、最后的input不能有./等符号的出现（不知为啥），也不能使用软连接，因此在这里选择直接拷贝过来。
#5、如果是双端测序，需要添加-p参数！！！
#6、-t 表示提取的type信息，如mRNA/exon/gene/ncRNA等，-g表示提取的ID类型，如gene_id,Name,gene等，但一定注意，两者要对上，即要保证特定type的annotation中有这类ID！！！
#7、注意这里的gff一定要与mapping时候的fa文件出处相同，不然会导致sort.bam匹配不上gff文件的基因位置！
