一、准备待分析文件
样本简况：
两个来自于化脓性链球菌的基因表达样本，每个样本有两个成对fastq文件,分别为 Read1 (R1) 和 Read2 (R2)。
样本一：(wild-type, 野生型)
RpSp_SP1_S5_L001_R1_001.2M.fastq
RpSp_SP1_S5_L001_R2_001.2M.fastq
样本二：(mutant,突变型)
RpSn_SP_S7_L001_R1_001.2M.fastq
RpSn_SP_S7_L001_R2_001.2M.fastq
（1）新建一个文件夹
mkdir RNAseq_bacterial
cd RNAseq_bacterial
将待比对文件拷贝进来。
（2）质量控制
fastp -i RpSp_SP1_S5_L001_R1_001.2M.fastq -I RpSp_SP1_S5_L001_R2_001.2M.fastq -o RpSp_SP1_S5_L001_R1_001.clean.fastq -O RpSp_SP1_S5_L001_R2_001.clean.fastq --trim_front1=5 --trim_front2=5 --trim_tail1=1 --trim_tail2=1 --html quality_report.html -w 1
fastp -i RpSn_SP_S7_L001_R1_001.2M.fastq -I RpSn_SP_S7_L001_R2_001.2M.fastq -o RpSn_SP_S7_L001_R1_001.clean.fastq -O RpSn_SP_S7_L001_R2_001.clean.fastq --trim_front1=5 --trim_front2=5 --trim_tail1=1 --trim_tail2=1 --html quality_report2.html -w 1
命令解读：
1.for single end data(not compressed):fastp -i in.fq -o out.fq
2.for paired end data (gzip compressed):fastp -i in.R1.fq.gz -I in.R2.fq.gz -o out.R1.fq.gz -O out.R2.fq.gz
3.the HTML report:-h或全称--html
4.the JSON report:-j
#Fastp支持全局修剪，这意味着修剪前面或尾部的所有读数。这个函数很有用，因为有时您想删除排序运行的一些周期。
5.--trim_front1和--trim_tail1:read1 or SE data开头和尾巴修剪，这里开头修剪五个碱基，尾巴修剪1个碱基。
6.--trim_front2和--trim_tail2:read2 of PE data开头和尾巴修剪，这里开头修剪五个碱基，尾巴修剪1个碱基。
#但是如果没有指定这些选项，它们将与read1选项相同，意味着trim_front2 = trim_front1和 trim_tail2 = trim_tail1。
7.-w:设置线程数，这里设置为1。
可参考fastp说明文件 https://github.com/OpenGene/fastp
（3）去除重复性的冗余序列
remove_redunt_pair_fsq.py RpSp_SP1_S5_L001_R1_001.clean.fastq RpSp_SP1_S5_L001_R2_001.clean.fastq 1 1unique RpSp_SP1_S5_L001_R1_001.clean.uniq.fastq RpSp_SP1_S5_L001_R2_001.clean.uniq.fastq
remove_redunt_pair_fsq.py RpSn_SP_S7_L001_R1_001.clean.fastq RpSn_SP_S7_L001_R2_001.clean.fastq 1 1unique RpSn_SP_S7_L001_R1_001.clean.uniq.fastq RpSn_SP_S7_L001_R2_001.clean.uniq.fastq
去除再不会用的文件。
rm RpSp_SP1_S5_L001_R1_001.clean.fastq
rm RpSp_SP1_S5_L001_R2_001.clean.fastq
rm RpSn_SP_S7_L001_R1_001.clean.fastq
rm RpSn_SP_S7_L001_R2_001.clean.fastq
（4）准备参考基因组文件
已准备好，且已格式化。
（5）利用比对软件 bwa 将短序列比对到参考基因组
bwa mem -t 2 -M AP53.fa RpSp_SP1_S5_L001_R1_001.clean.uniq.fastq RpSp_SP1_S5_L001_R2_001.clean.uniq.fastq > RpSp_SP1.sam
bwa mem -t 2 -M AP53.fa RpSn_SP_S7_L001_R1_001.clean.uniq.fastq RpSn_SP_S7_L001_R2_001.clean.uniq.fastq > RpSn_SP.sam
bwa men:使用的是bwa men算法。
-t:设置线程数。选择多线程的原因就是一个快字，使用多线程就是在正确的场景下通过设置正确个数的线程来最大化程序的运行速度。
-M:将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件。
参考：https://zhuanlan.zhihu.com/p/342250758
可参考bwa说明文档 http://bio-bwa.sourceforge.net/bwa.shtml
（6）利用 samtools 对 sam 格式的比对文件进行处理，以便进行后续分析
可参考samtools说明文档 http://www.htslib.org/doc/samtools-1.2.html
1.将 sam 格式转化成二进制格式 bam
samtools view -bt AP53.fa.fai -@ 2 -o RpSp_SP1.bam RpSp_SP1.sam 2>>samtools.log
samtools view -bt AP53.fa.fai -@ 2 -o RpSn_SP.bam RpSn_SP.sam 2>>samtools.log
命令解读：
-b:以bam文件格式输出，标准是以sam格式输出。
-t:使用一个list文件作为header的输入。
-T:使用序列fasta文件作为header的输入。
#-t和-T提供额外的参考资料。当sam输入文件没有@SQ headers时，这两个命令必须要用一个。
-@:添加额外的线程。
-o:指明输出文件名称。
命令 2>> 文件:将命令执行的错误输出结果重定向到指定的文件中，如果该文件中已包含数据，新数据将写入到原有内容的后面。这里就可以方便查看执行命令时是否发生错误和发生哪些错误。
2.对 bam 文件中的内容进行排序
samtools sort -O bam -@ 2 -o RpSp_SP1.sort.bam -T tmp_samtools RpSp_SP1.bam 2>>samtools.log
samtools sort -O bam -@ 2 -o RpSn_SP.sort.bam -T tmp_samtools RpSn_SP.bam 2>>samtools.log
命令解读：
-O:将文件以sam,bam或cram格式输出。默认情况下，samtools尝试基于-o文件名扩展名选择格式;如果输出是标准输出或不能推导出任何格式，则必须使用-O。这里指定以bam格式输出。
-@:设置排序和压缩线程的数量。默认情况下，操作是单线程的。
-o:指明输出文件名称。
-T:将临时文件写入PREFIX.nnnn.bam。这个选项是必需的。
命令 2>> 文件:将命令执行的错误输出结果重定向到指定的文件中，如果该文件中已包含数据，新数据将写入到原有内容的后面。这里就可以方便查看执行命令时是否发生错误和发生哪些错误。
3.为排序后的bam文件建立索引
samtools index RpSp_SP1.sort.bam
samtools index RpSn_SP.sort.bam
生成bai格式文件。
4.利用samtools对排序后的bam文件进行可视化
samtools tview RpSp_SP1.sort.bam AP53.fa
samtools tview RpSn_SP.sort.bam AP53.fa
samtools tview [-p chr:pos] [-s STR] [-d display] <in.sorted.bam> [ref.fasta]
5.利用IGV对排序后的bam文件进行可视化
（7）利用 htsep-count 计算比对到每个基因的短序列数目,其中参考基因组上的基因注释文件为gff格式
可参考 htseq-count 说明文档 ：https://htseq.readthedocs.io/en/release_0.11.1/count.html
参考基因组上的基因注释文件：AP53.rast.CDS.gff
htseq-count -f bam -s no -t CDS -i ID -m union --nonunique=none --secondary-alignments=ignore RpSp_SP1.sort.bam AP53.rast.CDS.gff > RpSp_SP1_expression.counts
htseq-count -f bam -s no -t CDS -i ID -m union --nonunique=none --secondary-alignments=ignore RpSn_SP.sort.bam AP53.rast.CDS.gff > RpSn_SP_expression.counts
命令解读：
-f:输入文件的格式。默认是sam格式。
-s:数据是否来自a strand-specific assay。如果是no，一个read 被认为与一个特征重叠，不管它是映射到与该特征相同还是相反的链上。如果是yes，单端读取，读取必须映射到与特征相同的链上；双端读取，第一次读取必须在同一条链上，第二次读取在相反的链上。如果是reverse，这些读取的规则相反。
-t:使用的特征类型(GFF文件的第三列)，其他类型的所有特征都被忽略(默认，适用于使用ensemble GTF文件的RNA-Seq分析:外显子)。
-i:GFF属性特征ID。具有相同特征ID的多个GFF线将被视为同一特征的一部分。特性ID用于标识输出表中的计数。默认值是gene_id，适合使用ensemble GTF文件进行RNA-Seq分析。
-m:模式处理重叠多个特性的读取。模式是联合、相交-严格和相交-非空(默认为union联合)。
--nonunique:模式来处理与重叠模式中的多个特性对齐或分配给该特性的读取。nonunique是none和all(默认值:none)。
--secondary-alignments:处理辅助对齐的模式(SAM标志0x100)。可以是score和ignore(默认值:score)。
（8）采用RPKM方法生成归一化因子，以便进行表达水平归一化
python3 /bioinfo/script/rpkm_normalization_factor.py -counts AP53_counts.txt -gff AP53.rast.gff -out AP53_rpkmFactor.txt
（9）基因差异表达计算
可参考说明文件：https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
1.执行命令R 进入R环境,并读取差异表达分析包 DESeq2
R
library(DESeq2)
2.读取短片段比对的基因计数文件 AP53_counts.txt 和归一化因子文件 AP53_rpkmFactor.txt,并查看其内容
cuntTab2 <- read.delim("AP53_counts.txt",header=TRUE,row.names=1,stringsAsFactors=TRUE)
normFactors <- data.matrix(read.table("AP53_rpkmFactor.txt",header=TRUE,row.names=1))
head(cuntTab2)
head(normFactors)
3.样品条件设置以及差异分析所需数据框生成
conds2 <- c("RpSp","RpSn","RpSp","RpSn")
samples <- data.frame(row.names=colnames(cuntTab2),conds2)
normTab <- DESeqDataSetFromMatrix(countData=cuntTab2,colData=samples,design=~conds2)
length(normTab)
4.基因表达计数的归一化 normalization,并作boxplot图
normalizationFactors(normTab) <- normFactors
plotCounts(normTab,'AP53_1668',intgroup='conds2')
normCunt <- counts(normTab,normalized=TRUE)
boxplot(normCunt,ylim=c(-1,2000))
5.差异表达计算,并写入文件
normDiff <- DESeq(normTab)
resDiff <- results(normDiff)
head(resDiff)
write.table(resDiff,"AP53_DESeq2_withNorm.txt",row.names=TRUE,quote=FALSE,sep="\t")
6.作图查看结果: p-value分布图, MA图, 火山图
hist(resDiff$pvalue,breaks=10,col="grey",xlab="p-value")
plot(resDiff$baseMean,resDiff$log2FoldChange,log='x',col=ifelse(resDiff$baseMean >= 100 & abs(resDiff$log2FoldChange) >= 2,"red","black"),xlab="mean normalized expression",ylab="log2FoldChange")
plot(resDiff$baseMean,resDiff$padj,log='x',col=ifelse(resDiff$baseMean >= 100 & abs(resDiff$padj) <= 0.05,"red","black"),xlab="mean normalized expression",ylab="p-vlaue")
plot(resDiff$log2FoldChange,-log2(resDiff$padj),col=ifelse(abs(resDiff$log2FoldChange) >= 2 & abs(resDiff$padj) <= 0.05,"red","black"),xlab="log2FoldChange",ylab="-log2Pvalue")
