1. 软件安装
1.1需要更改镜像源配置
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
1.2 安装软件
创建环境
conda create -n rnaseq python=3 #创建名为rna的软件安装环境
conda info --envs #查看当前conda环境
source activate rnaseq #激活conda的rna环境
转录组分析需要用到的软件列表
质控
fastqc , multiqc, trimmomatic, cutadapt, trim-galore
比对
star, hisat2, bowtie2, tophat, bwa, subread
计数
htseq, bedtools, deeptools, salmon
安装软件
conda install -y sra-tools
conda install -y trimmomatic
conda install -y cutadapt multiqc
conda install -y trim-galore
conda install -y star hisat2 bowtie2
conda install -y subread tophat htseq bedtools deeptools
conda install -y salmon
source deactivate #注销当前的rna环境
2 数据下载
2.1.1 下载SRA数据
从GEO数据库下载到SRR_Acc_List.txt的文档，文档内容是SRR号码（一个号码占据一行）
#创建目录
mkdir -p project/{sra,fastq,rmdup,clean,counts,align,peaks,motif,qc/{raw,trimed}}
cd project/
#或者自己制作list文件，将下列数据写入vim文件中
vim SRR_Acc_List.txt
SRR1266976
SRR1266977
SRR1266978
SRR1266979
SRR1266980
SRR1266981
cat SRR_Acc_List.txt | while read id; do (prefetch ${id} );done
2.1.2 Aspera Connect下载
SRA数据库下载：数据的存放地址是
ftp://ftp.sra.ebi.ac.uk/vol1
举例：下载SRR1577019文件，首先我需要找到地址，去ftp://ftp.sra.ebi.ac.uk/vol1，一层层寻找，直至找到，ftp://ftp.sra.ebi.ac.uk/vol1/srr/SRR157/009/SRR1577019
去geo官网搜索srr然后下载获取SRR_Acc_List.txt 然后观察是err+7位数的，所以直接使用7位数的代码。（7位数y取最后一位，8位数y取最后2位，6位数把y变量删掉，应该是这样。）一般来说，NCBI的sra文件前面的地址都是一样的 ~/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk，那么写脚本批量下载也就不难了！最后那个 ./ 表示下载后的路径。
- 1批量下载srr文件（推荐）
cd ~/
mkdir -p rnaseq/{sra,fastq,rmdup,clean,align,peaks,qc/{raw,trimed}}
cd ~/rnaseq/
#创建文件夹
source activate rnaseq
cat SraAccList.txt | while read id
do
x=$(echo $id | cut -b1-6)
y=$(echo $id | cut -b10-10)
echo $id
ascp -QT -k 1 -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/srr/$x/00$y/$id/ ~/rnaseq/sra/
done
下载单个srr文件
ascp -QT -k 1 -l 500m -P 33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/srr/SRR126/000/SRR1266980 ./sra
下载fastq文件
cat SRR_Acc_List.txt | while read id
do
x=$(echo $id | cut -b1-6)
y=$(echo $id | cut -b10-10)
echo $id
ascp -QT -k 1 -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/$x/00$y/$id/ ./sra
done
如果下载的fastq数据是2个文件 *_1.fastq 和 *_2.fastq表明是双端测序。
Aspera用法如下：
ascp [参数] 目标文件 保存路径
-v verbose mode 实时知道程序在干啥
-T 取消加密，否则有时候数据下载不了
-i 提供私钥文件的地址
-l 设置最大传输速度，一般200m到500m，如果不设置，反而速度会比较低，可能有个较低的默认值
-k 断点续传，一般设置为值1
-Q 一般加上它
-P 提供SSH port，端口一般是33001
3. sra文件转fastq文件：
- sra转化为fastq文件可以使用sratoolkit中的fastq-dump命令。
fastq-dump --split-3 filename
其中 --split-3 参数代表着如果是单端测序就生成一个 、.fastq文件，如果是双端测序就生成_1.fastq 和_2.fastq 文件；即单端测序输出1个fastq文件，双端测序输出2个fastq文件。
进入到sra文件中,将SRR文件夹去掉，数据文件直接保存在sra目录中，我们可以用下述代码进行批量的格式转化：
首先制作如下对应的文本
SRR8444250 ES_Input
SRR8444251 ES_Smad2_ChIPseq
SRR8444256 EB_AC_Input
SRR8444257 EB-SB_Smad2_ChIPseq
SRR8444258 EB-AC_Smad2_ChIPseq
然后
## 下面需要用循环
source activate rnaseq
cd ~/rnaseq/sra/
#把list.txt文件复制到sra目录下
## 下面用到的 list.txt 文件，就是上面自行制作的对应的文件。
cat SraAccList.txt | while read id
do
echo $id
arr=($id)
srr=${arr[0]} #通过中括号获取数组中的某一个元素：${arr[0]} 得到第一个元素，${arr[0]} 第二个
sample=${arr[1]}
(nohup fastq-dump -A $sample -O ~/rnaseq/fastq/ --gzip --split-3 $srr &)
done
#-A -$sample是输出文件名字
#--gzip打包节省空间，且对后续分析不影响，如果有影响就把这个参数去掉，-O输出到目录
4. 检测数据质量及质量过滤
4.1.1 检测数据质量
fastqc生成质控报告
cd ~/rnaseq/qc/raw/
files=~/rnaseq/fastq/*.fastq.gz
ls $files | while read id
do
echo $id
fastqc $id -t 10 -o ~/rnaseq/qc/raw/
done
#-t：调用核心数目
#-q：安静运行，运行过程中不会生成报告，在结束时将报告生成一个文件
#-o ./：文件输出当前目录
4.1.2 multiqc将各个样本的质控报告整合为一个
cd ~/rnaseq/qc/raw/
files=~/rnaseq/qc/raw/*.zip
multiqc *.zip
4.2 使用trim_galore软件对fastq文件进行质量控制-过滤
4.2.1 首先判断phred 33 和phred 64
如何判断phred 33 和phred 64?
有时候得到的原始fastq文件，无法知道质量值体系，你就无法进行质量值的过滤，我们可以在正常情况下，按照上面的表对回去，统计一下几条reads的最大和最小质量值的区间，就可以知道到底是phred 33 还是phred 64体系。
根据上图中Phred+33与Phred+64所使用的质量字符范围的不同，可以对fastq文件中质量得分的编码方式做出判断(第四行是测序质量，查看第四行符号是在Phred+64还是Phred+33范围)。图中显示，ASCII值小于等于58（相应的质量得分小于等于25）对应的字符只有在Phred+33的编码中被使用，所有Phred+64所使用的字符的ASCII值都大于等于59。在通常情况下，ASCII值大于等于74的字符只出现在Phred+64中。利用这些信息即可在程序中进行判断。（近几年应该都是Phred+33）
4.2.2 trim_galore处理fastq数据
批量处理双端测序数据
cd ~/rnaseq/clean/
files=~/rnaseq/fastq/*.gz
ls $files | while read id
do
sample=$(basename $id)
fq=${sample%%_?.fastq*}
echo $fq
trim_galore -q 10 --phred33 --length 25 --stringency 4 --paired ~/rnaseq/fastq/${fq}_1.fastq.gz ~/rnaseq/fastq/${fq}_2.fastq.gz --gzip -o ~/rnaseq/clean/
done
批量处理单端测序数据
cd ~/rnaseq/clean
files=~/rnaseq/fastq/*.gz
ls $files | while read id
do
sample=$(basename $id)
echo $sample
trim_galore -q 10 --phred33 --length 25 --stringency 4 ~/rnaseq/fastq/${sample} -o ~/rnaseq/clean/
done
#具体稍作调整
- fastqc对质控后的数据检测
cd ~/rnaseq/qc #切换到qc目录
#比较经trim_galore处理前后数据质量
ls ~/rnaseq/clean/*.gz | xargs fastqc -t 10 -o ~/rnaseq/qc/trimed/
#-t：调用核心数目
#-q：安静运行，运行过程中不会生成报告，在结束时将报告生成一个文件
#-o ./：文件输出当前目录
#*.gz：输入文件
- multiqc可以对几个fastqc报告文件进行总结并汇总到一个报告文件中，以更直观到防止展示。使用方法
multiqc <analysis directory>
multiqc ~/rnaseq/qc/trimed/ -o ~/rnaseq/qc/trimed/
# multiqc 后面直接跟上 fastqc 结果所在的文件夹
# -n 输出结果的文件名前缀
# -o 设置输出结果的文件夹
#(有时间再看)
5. hisat2 比对
官方手册：https://daehwankimlab.github.io/hisat2/manual/
参考基因组下载
有三大全文网站提供参考基因组下载，它们分别是：
1.NCBI （https://www.ncbi.nlm.nih.gov/grc）
2.UCSC (http://hgdownload.soe.ucsc.edu/downloads.html)
3.Ensemble （http://asia.ensembl.org/index.html?redirect=no）
目前最常用的人和小鼠的参考基因组版本如下（Jimmy总结）：
|NCBI | UCSC| Ensemble|
|GRCh36 | hg18 | ENSEMBL release_52 |
|GRCh37 | hg19 | ENSEMBL release_59/61/64/68/69/75|
|GRCh38 | hg38 | ENSEMBL release_76/77/78/80/81/82|
5.1 建立index
5.1.1自己建立index
#下载基因组文件
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz（UCSC版本基因组）
(我的已下载，存储在~/biosoft/Reference genome/UCSC/目录下是“chromFa.tar.gz”)
# 解压缩
tar -zxvf chromFa.tar.gz
#解压出来的新文件都是同一类型：chr*.fa。
chr1.fa
chr10.fa
chr11.fa
chr11_gl000202_random.fa
chr12.fa
chr13.fa
[...]
cat *.fa > hg19.fa #将解压出来的所有.fa结尾的文件合并为一个文件，即hg19.fa，到此为止，参考基因组就准备好了。
hisat2-build -p 4 hg19.fa genome #建立HISAT2索引文件
#下载基因组注释文件
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/genes/hg19.ncbiRefSeq.gtf.gz
gunzip hg19.ncbiRefSeq.gtf.gz
5.1.2hisat 官网下载index
参考
index文件直接去官网下载
(http://daehwankimlab.github.io/hisat2/download/)，如下，为下载的小鼠的index。视频教程
解压文件，解压过程中会在当前文件夹下创建mm10文件，解压后的文件就在mm10文件夹中。
tar -zxvf /data/mouse_genome/mm10_genome.tar.gz
#我的已经下载，存储在~/biosoft/reference_hisat2/index/目录下
5.2 hisat2 比对
双端测序
cd ~/rnaseq/align
#定义索引路径
index=~/biosoft/Reference_genome/hisat2_index/ensembl/mm10/genome
#一定要注意 ~/biosoft/reference_hisat2/index/mm10/genome文件的目录，genome这个不是文件夹，是index文件的前缀
files=~/rnaseq/clean/*_1.fq.gz
outputdir=~/rnaseq/align/
ls $files | while read id
do
sample=${id%%_?_val_?.fq.gz}
files_name=$(basename $sample)
echo ${files_name}
hisat2 -p 10 --dta -x $index -1 ${sample}_1_val_1.fq.gz -2 ${sample}_2_val_2.fq.gz | samtools sort -@ 5 -o $outputdir/${files_name}.sorted.bam && samtools index $outputdir/${files_name}.sorted.bam $outputdir/${files_name}.sorted.bam.bai && samtools flagstat $outputdir/${files_name}.sorted.bam >${files_name}.flagstat.txt
done
或者用下述代码
fq1=~/rnaseq/clean/*_1.fq.gz
fq2=~/rnaseq/clean/*_2.fq.gz
outputdir=~/rnaseq/align/
ls $files_1 | while read id
do
fq1=$id
for i in `ls $files_2 | grep "${fq1%%_1_val*}"`
do
fq2=$i
file=$(basename $fq1)
fq=${file%%_1_val*}
done
echo $fq1
echo $fq2
hisat2 -p 10 --dta -x $index -1 $fq1 -2 $fq2 | samtools sort -@ 5 -o $outputdir/${fq}.sorted.bam && samtools index $outputdir/${fq}.sorted.bam $outputdir/${fq}.sorted.bam.bai && samtools flagstat $outputdir/${fq}.sorted.bam >$fq.flagstat.txt
done
#-x 参考基因组索引文件的前缀。
#-1 双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
#-2 双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
#-U 单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。
#–sra-acc 输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。
#-S 指定输出的SAM文件。
#-t/--time 输出搜索阶段所花费的wall-clock时间 print wall-clock time taken by search phases
#另外一定要加上 --dta，后续用Stringtie处理会更容易一些，这StringTie
#使用说明里面的一句话：【NOTE: be sure to run HISAT2 with the --dta
#option for alignment, or your results will suffer.】
# 常用参数
hisat2 [option] -1 <m1> -2 <m2> -x <index_base> -S <out.sam>
-x <path/to/base> #索引的前缀,也可以添加路径，例如 ～/genome/homo_GRCh38.dna.primary
-1 #双端测序的第一个文件，若有多组，则用逗号隔开，名字与2要匹配，如-1 flyA_1.fq,flyB_1.fq
-2 #类比上，名字与1要匹配，如-2 flyA_2.fq,flyB_2.fq
-U #单端数据文件
-S #保存到的sam文件，默认输出到标准输出
# 以下为参数
-p #线程数
--dta # > 在下游使用stringtie组装的时候量身定做的参数。使用此选项，
# ^ HISAT2需要更长的锚长度才能从头发现拼接位点，这样可以减少与短锚的对齐，
# ^ 这有助于转录汇编程序显着提高计算和内存使用率。
# ^ 当然下游不使用stringtie也可以使用此参数减少短锚定read
--dta-cufflinks #下游使用cufflinks则需要添加此参数
-q #指定读取的测序文件是fastq文件（含有质量信息）
-f #指定读取的测序文件是fa文件
-t #打印加载索引文件和对齐读数所需的时钟时间
--phred33 #质量表示方法，如果使用fq文件，建议选择，同理有小众的--phred64
--min-intronlen #设置最小内含子的长度，默认值20
--max-intronlen #设置最大内含子长度，默认500000
--known-splicesite-infile <path> #提供已知的剪接位点列表，HISAT2利用这些列表比对,可以用hisat2_extract_splice_sites.py和gtf生成信息
--novel-splicesite-outfile <path> #在结果中生成一个剪接位点的列表
--novel-splicesite-infile <path> # > 可以使用novel-splicesite-outfile所生成的剪接列表作为
# ^ 新剪接列表，应该可以自己定义。为特定基因。
单端
source activate rnaseq
cd ~/rnaseq/align/
index=~/biosoft/reference_hisat2/index/mm10/genome
inputdir=~/rnaseq/clean/*fq.gz
outputdir=~/rnaseq/align/
ls $inputdir | while read id
do
sample=$(basename -s _trimmed.fq.gz $id)
echo $sample
hisat2 -p 10 -x $index -U $id | samtools sort -@ 5 -o $outputdir/${sample}.sorted.bam && samtools index $outputdir/${sample}.sorted.bam $outputdir/${sample}.sorted.bam.bai && samtools flagstat $outputdir/${sample}.sorted.bam >$sample.flagstat.txt
done
sam转bam文件
samtools主要参数：
-view BAM-SAM/SAM-BAM转换和提取部分比对
-sort 比对排序
-merge 聚合多个排序比对
-index 索引排序比对
-faidx 建立FASTA索引，提取部分序列
5.3 加载IGV，可视化比对结果
载入参考序列，注释和BAM文件?学IGV必看的初级教程
samtools 软件进行格式转换
SAM文件和BAM文件 samtools 是针对比对回帖的结果——sam和bam格式文件的进一步分析使用的软件。sam格式文件由于体量过大，一般都是使用bam文件来进行存储。由于bam文件是二进制存储所以文件大小比sam格式文件小许多，大约是sam格式体积的1/6 。 samtools将sam转换bam文件
samtools view -S seq.sam -b > seq.bam #文件格式转换
samtools sort seq.bam -0 seq_sorted.bam ##将bam文件排序
samtools index seq_sorted.bam #对排序后对bam文件索引生成bai格式文件，用于快速随机处理。
至此一个回帖到基因组对RNA-seq文件构建完成。这个seq_sourted.bam文件可以通过samtools或者IGV( Integrative Genomics Viewer)独立软件进行查看。在IGV软件中载入seq_sourted.bam文件。 可以很直观清晰地观察到reads在基因组中的回帖情况和外显子与内含子的关系。
6 计算count
计算RNA-seq测序reads对在基因组中对比对深度。 计数工具：feature counts 公式构建：
feature counts -T 6 -t exon -g gene_id -a <gencode.gtf> -o seq_featurecount.txt <seq.bam>
用户手册下载：http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf
featurecounts说明在P31
SAM/BAM文件和GTF/GFF/SAF文件需要来自同一个参考基因组，即必须参考基因组和GTF/GFF/SAF文件来自同一个网站，同一个版本
- 参数：
-a 输入GTF/GFF基因组注释文件
-p 这个参数是针对paired-end数据
-F 指定-a注释文件的格式，默认是GTF
-g 从注释文件中提取Meta-features信息用于read count，默认是gene_id
-t 跟-g一样的意思，其是默认将exon作为一个feature
-o 输出文件
-T 多线程数
(gencode注释文件对应版本说明说明)https://www.jianshu.com/p/eda5263d4494
gtf下载（ensembl版本 地址）
gtf下载gencode版
gtf下载（UCSC版本 地址）https://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/genes/mm10.refGene.gtf.gz（小鼠mm10）
如下图：https://hgdownload.soe.ucsc.edu/进入UCSC主页面后点download按钮按照下图进行（refGene.gtf.gz就是UCSC版本的gtf文件）
参考：https://www.bilibili.com/read/cv10447213/
source activate rnaseq
cd ~/rnaseq/counts/
#需要下载gtf文件（我用ensembl版）
gtf=~/biosoft/Reference_genome/gtf/ensembl/mm10/Mus_musculus.GRCm38.102.gtf.gz
inputdir=~/rnaseq/align/*.bam
featureCounts -T 10 -t exon -g gene_id -a $gtf $inputdir -o all.id.txt
#-g # 注释文件中提取对Meta-feature 默认是gene_id, -t # 提取注释文件中的Meta-feature
#默认是 exon, -p #参数是针对paired-end 数据, -a #输入GTF/GFF 注释文件 -o #输出文件
# 对定量结果质控
multiqc all.id.txt.summary
# 得到表达矩阵
cat all.id.txt | cut -f1,7- > counts.txt
source deactivate
featureCounts [options] <input.file>
<input.file> #输入文件，sam/bam 支持多个文件输入
-a <gtf> #参考gtf注释文件，支持gz压缩
-F <参考文件格式> #默认GTF
-T <int> #线程数 1-32
-J #对可变剪接计数
-G <str> #有-J设置时，用来提供参考基因组辅助寻找可变剪接
-M #如果设置了-M则多重比对会被统计
-o <out.file> #输出文件的名字，输出文件的内容是read的统计数目
-O #允许多重比对，当一个read比对到多个feature或多个metafeature的时候，这条read会被统计多次
-p #声明测序类型是paired-end，此时，后续统计fragment而不是read
-B #在-p设置时，只有两端都比对上的fragment才会被统计
-C #在-p设置时，-C声明比对到不同染色体的fragment不计数
-d <int> #最短的fragment ，默认50
-D <int> #最长的fragment，默认600
-f #设置后会统计feature水平数据，如exon-level，否则会统计meta-feature层面数据，如gene-level，（？应该是和-g冲突,一般与-t连用?）
-g <str> #参考的gtf提供的时候，我们需要提供一个id identifier 来将feature水平的统计汇总为meta-feature水平的统计，默认为gene_id，可以选择transcript_id、gene_name、transcript_name等
-t <str> #设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon”（？应该和-g冲突）
