a.整体流程及相关软件
数据资源下载，参考基因组及参考注释组
数据质控及相关软件
fastp (代替fastqc，multiqc， trimmomatic，cutadapt,trim_galore)比对及相关软件
HISAT2，STAR，TOPHAT2，bowtie2，bwa...
本文用HISAT2统计及相关软件
featureCounts，htseq-counts，bedtools...均一化及差异分析
DEseq2， edgeR，limma，cuffdiff...
b.环境搭建及软件安装
conda create -n rnaseq python=3.8 r-base=4.0
#搭建conda环境，指定python及R版本
conda activate rnaseq
#激活conda环境
conda install -c bioconda fastp
#conda install -c bioconda fastqc
#conda install -c bioconda multiqc
#conda install -c bioconda trimmomatic
#conda install -c bioconda cutadapt
#安装质控相关软件，#号为可选软件，推荐用fastp代替，方便
conda install -c bioconda hisat2
conda install -c bioconda samtools
conda install -c bioconda subread
conda install -c bioconda htseq
conda install -c bioconda bioconductor-deseq2
conda install -c bioconda bioconductor-edger
c.具体流程
1.数据质控和过滤
本文中使用fastp一步到位对数进行质控、过滤并统计。
具体参数见前文fastp部分。测序数据过滤原理相同（去接头、去低质量等）。
usage:fastp -i <in1> -o <out1> [-I <in1> -O <out2>] [options...]
#1 PE测序
fastp -i sample.R1.fastq.gz -o sample.R1.clean.fq.gz \
-I sample.R2.fastq.gz -O sample.R2.clean.fq.gz \
-5 --cut_front_window_size 4 --cut_front_mean_quality 20 \
-3 --cut_tail_window_size 4 --cut_tail_mean_quality 20 \
--cut_right --cut_right_window_size 4 --cut_right_mean_quality 20 \
--detect_adapter_for_pe \
-q 15 -u 40 -e 20 -n 5 -l 30 -p -P 20 -w 4 -R "sample"
#2 SE测序
fastp -i sample.fq -o sample.clean.fq \
-5 --cut_front_window_size 4 --cut_front_mean_quality 20 \
-3 --cut_tail_window_size 4 --cut_tail_mean_quality 20 \
--cut_right --cut_right_window_size 4 \
--cut_right_mean_quality 20 \
-f 5 -q 15 -u 40 -e 20 -n 5 -p -P 20 -w 4 -R "sample"
过滤完成后用fastqc再跑一遍看质量
fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
-o --outdir FastQC生成的报告文件的储存路径，生成的报告的文件名是根据输入来定的
-f --format 指定输入文件的格式
--extract 生成的报告默认会打包成1个压缩文件，使用这个参数是让程序不打包
-t --threads 选择程序运行的线程数，每个线程会占用250MB内存，越多越快咯
--min_length 设置序列的最小长度，≥最长read的长度
-c --contaminants 污染物选项，输入的是一个文件，格式是Name [Tab] Sequence，里面是可能的污染序列，如果有这个选项，FastQC会在计算时候评估污染的情况，并在统计的时候进行分析，一般用不到
-a --adapters 也是输入一个文件，文件的格式Name [Tab] Sequence，储存的是测序的adpater序列信息，如果不输入，目前版本的FastQC就按照通用引物来评估序列时候有adapter的残留
-q --quiet 安静运行模式，一般不选这个选项的时候，程序会实时报告运行的状况
跑完fastqc后用multiqc汇总质量报告，网页查看
multiqc ./
2 HISAT2比对
##构建索引
hisat2-build -p 20 ~/ref/hg19.fa ./hg19/hg19
#-p 线程数
#~/ref/hg19.fa 参考基因组
#./hg19/hg19 输出文件及前缀
##比对
hisat2 -x ~/ref/hisat2_index/hg19/hg19 -p 40 --dta -1 ../2.cleandata/$i.1.clean.fq.gz -2 ../2.cleandata/$i.2.clean.fq.gz -S $i.sam
#-x 构建好的参考基因组索引
#-p 线程数
#--dta 用于组装需要的参数
#-1 第一条reads
#-2 第二条reads
#-S 输出为sam文件
##sam批量转换为bam
ls *.sam |while read id;do (samtools view -@ 5 -b -S ${id} -o $(basename ${id} ".sam").bam);done
##bam文件排序
samtools sort -n sample.bam -o sample.sort.bam
##bam文件批量质控
ls *.bam |while read id;do (samtools flagstat -@ 5 ${id} > $(basename ${id} ".bam").flagstat);done
##Featurecounts做成表达矩阵（单个样品）
featureCounts -T 5 -p -t exon -g gene_id -a ~/ref/hg19.ensGene.gtf -o $i.enscounts.txt $i.sort.bam
##Featurecounts做成表达矩阵（所有样品）
featureCounts -T 5 -p -t exon -g gene_id -a ~/ref/hg19.ensGene.gtf -o all.enscounts.txt *.sort.bam
#说明：很奇怪，用refGene.gtf 每个样本能找出28267条信息，但是用ensGene.gtf每个样本能找出60236条记录，不知道发生了什么，可能是ensemble数据库基因注释较多。我后续用的ens注释做分析。
#详细参数
#input file 输入的bam/sam文件，支持多个文件输入
-a < string > 参考gtf文件名，支持gz压缩文件
-A 提供一个逗号分割为两列的文件，一列为gtf中的染色体名，另一列为read中对应的染色体名，用于将gtf和read中的名称进行统一匹配，注意该文件提交时不需要列名
-J 对可变剪切进行计数
-G < string > 当-J设置的时候，通过-G提供一个比对的时候使用的参考基因组文件，辅助寻找可变剪切
-M 如果设置-M，多重map的read将会被统计到
-o < string > 输出文件的名字，输出文件的内容为read 的统计数目
-O 允许多重比对，即当一个read比对到多个feature或多个metafeature的时候，这条read会被统计多次
-T 线程数目，1~32
下面是有关featrue/metafeature选择的参数 参数说明
-p 只能用在paired-end的情况中，会统计fragment而不统计read
-B 在-p选择的条件下，只有两端read都比对上的fragment才会被统计
-C 如果-C被设置，那融合的fragment（比对到不同染色体上的fragment）就不会被计数，这个只有在-p被设置的条件下使用
-d < int > 最短的fragment，默认是50
-D < int > 最长的fragmen，默认是600
-f 如果-f被设置，那将会统计feature层面的数据，如exon-level，否则会统计meta-feature层面的数据，如gene-levels
-g < string > 当参考的gtf提供的时候，我们需要提供一个id identifier 来将feature水平的统计汇总为meta-feature水平的统计，默认为gene_id，注意！选择gtf中提供的id identifier！！！
-t < string > 设置feature-type，-t指定的必须是gtf中有的feature，同时read只有落到这些feature上才会被统计到，默认是“exon”
##Featurecounts结果说明
Geneid Chr Start End Strand Length ck1.sort.bam
#第一列 Geneid
#第二列 chr，可能对应多个位置
#第三列 起始位点，可对应多个位点
#第四列 终止位点，跟前边一一对应
#第五列 +-链信息
#第六列 基因长度
#第七列 样本表达信息（最重要的一列）
##提取信息做成数据表达矩阵
for i in $(cat ../sample.txt);do awk -F "\t" '{print $1"\t"$7}' >$i.expr.txt $i.enscounts.txt;done
#提取第1列（geneid）和第7列（expression）存入$i.expr.txt
ls *.expr.txt|while read dd;do sed -i '1d' $dd;done
#将第一行信息删除
剩下的用R语言统计比较方便了
3 差异表达基因统计
说明文件请参考RNA-seq(6): reads计数，合并矩阵并进行注释
计数分为三个水平： gene-level, transcript-level, exon-usage-level
标准化方法： FPKM RPKM TMM TPM
该文使用DEseq2
setwd("./")
