a.整體流程及相關(guān)軟件

1. 數(shù)據(jù)資源下載，參考基因組及參考注釋組

2. 數(shù)據(jù)質(zhì)控及相關(guān)軟件
   fastp (代替fastqc，multiqc， trimmomatic，cutadapt,trim_galore)

3. 比對(duì)及相關(guān)軟件
   HISAT2，STAR，TOPHAT2，bowtie2，bwa...
   本文用HISAT2

4. 統(tǒng)計(jì)及相關(guān)軟件
   featureCounts，htseq-counts，bedtools...

5. 均一化及差異分析
   DEseq2， edgeR，limma，cuffdiff...

b.環(huán)境搭建及軟件安裝

conda create -n rnaseq python=3.8 r-base=4.0
#搭建conda環(huán)境，指定python及R版本

conda activate rnaseq
#激活conda環(huán)境

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
#安裝質(zhì)控相關(guān)軟件，#號(hào)為可選軟件，推薦用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.數(shù)據(jù)質(zhì)控和過濾

本文中使用fastp一步到位對(duì)數(shù)進(jìn)行質(zhì)控、過濾并統(tǒng)計(jì)。
具體參數(shù)見前文fastp部分。測序數(shù)據(jù)過濾原理相同（去接頭、去低質(zhì)量等）。

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再跑一遍看質(zhì)量

fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN

-o --outdir       FastQC生成的報(bào)告文件的儲(chǔ)存路徑，生成的報(bào)告的文件名是根據(jù)輸入來定的
-f --format       指定輸入文件的格式
--extract         生成的報(bào)告默認(rèn)會(huì)打包成1個(gè)壓縮文件，使用這個(gè)參數(shù)是讓程序不打包
-t --threads      選擇程序運(yùn)行的線程數(shù)，每個(gè)線程會(huì)占用250MB內(nèi)存，越多越快咯
--min_length      設(shè)置序列的最小長度，≥最長read的長度
-c --contaminants 污染物選項(xiàng)，輸入的是一個(gè)文件，格式是Name [Tab] Sequence，里面是可能的污染序列，如果有這個(gè)選項(xiàng)，F(xiàn)astQC會(huì)在計(jì)算時(shí)候評(píng)估污染的情況，并在統(tǒng)計(jì)的時(shí)候進(jìn)行分析，一般用不到
-a --adapters     也是輸入一個(gè)文件，文件的格式Name [Tab] Sequence，儲(chǔ)存的是測序的adpater序列信息，如果不輸入，目前版本的FastQC就按照通用引物來評(píng)估序列時(shí)候有adapter的殘留
-q --quiet        安靜運(yùn)行模式，一般不選這個(gè)選項(xiàng)的時(shí)候，程序會(huì)實(shí)時(shí)報(bào)告運(yùn)行的狀況

跑完fastqc后用multiqc匯總質(zhì)量報(bào)告，網(wǎng)頁查看

multiqc ./

2 HISAT2比對(duì)

##構(gòu)建索引
hisat2-build -p 20 ~/ref/hg19.fa ./hg19/hg19 
#-p 線程數(shù)
#~/ref/hg19.fa 參考基因組
#./hg19/hg19 輸出文件及前綴

##比對(duì)
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 構(gòu)建好的參考基因組索引
#-p 線程數(shù)
#--dta 用于組裝需要的參數(shù)
#-1 第一條reads
#-2 第二條reads
#-S 輸出為sam文件

##sam批量轉(zhuǎn)換為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文件批量質(zhì)控
ls *.bam |while read id;do (samtools flagstat -@ 5  ${id} > $(basename ${id} ".bam").flagstat);done

##Featurecounts做成表達(dá)矩陣（單個(gè)樣品）
featureCounts -T 5 -p -t exon -g gene_id -a ~/ref/hg19.ensGene.gtf -o $i.enscounts.txt $i.sort.bam

##Featurecounts做成表達(dá)矩陣（所有樣品）
featureCounts -T 5 -p -t exon -g gene_id -a ~/ref/hg19.ensGene.gtf -o all.enscounts.txt *.sort.bam

#說明：很奇怪，用refGene.gtf 每個(gè)樣本能找出28267條信息，但是用ensGene.gtf每個(gè)樣本能找出60236條記錄，不知道發(fā)生了什么，可能是ensemble數(shù)據(jù)庫基因注釋較多。我后續(xù)用的ens注釋做分析。
#詳細(xì)參數(shù)
#input file 輸入的bam/sam文件，支持多個(gè)文件輸入
-a < string >   參考gtf文件名，支持gz壓縮文件
-A  提供一個(gè)逗號(hào)分割為兩列的文件，一列為gtf中的染色體名，另一列為read中對(duì)應(yīng)的染色體名，用于將gtf和read中的名稱進(jìn)行統(tǒng)一匹配，注意該文件提交時(shí)不需要列名
-J  對(duì)可變剪切進(jìn)行計(jì)數(shù)
-G < string >   當(dāng)-J設(shè)置的時(shí)候，通過-G提供一個(gè)比對(duì)的時(shí)候使用的參考基因組文件，輔助尋找可變剪切
-M  如果設(shè)置-M，多重map的read將會(huì)被統(tǒng)計(jì)到
-o < string >   輸出文件的名字，輸出文件的內(nèi)容為read 的統(tǒng)計(jì)數(shù)目
-O  允許多重比對(duì)，即當(dāng)一個(gè)read比對(duì)到多個(gè)feature或多個(gè)metafeature的時(shí)候，這條read會(huì)被統(tǒng)計(jì)多次
-T  線程數(shù)目，1~32
下面是有關(guān)featrue/metafeature選擇的參數(shù)   參數(shù)說明
-p  只能用在paired-end的情況中，會(huì)統(tǒng)計(jì)fragment而不統(tǒng)計(jì)read
-B  在-p選擇的條件下，只有兩端read都比對(duì)上的fragment才會(huì)被統(tǒng)計(jì)
-C  如果-C被設(shè)置，那融合的fragment（比對(duì)到不同染色體上的fragment）就不會(huì)被計(jì)數(shù)，這個(gè)只有在-p被設(shè)置的條件下使用
-d < int >  最短的fragment，默認(rèn)是50
-D < int >  最長的fragmen，默認(rèn)是600
-f  如果-f被設(shè)置，那將會(huì)統(tǒng)計(jì)feature層面的數(shù)據(jù)，如exon-level，否則會(huì)統(tǒng)計(jì)meta-feature層面的數(shù)據(jù)，如gene-levels
-g < string >   當(dāng)參考的gtf提供的時(shí)候，我們需要提供一個(gè)id identifier 來將feature水平的統(tǒng)計(jì)匯總為meta-feature水平的統(tǒng)計(jì)，默認(rèn)為gene_id，注意！選擇gtf中提供的id identifier！?。?-t < string >   設(shè)置feature-type，-t指定的必須是gtf中有的feature，同時(shí)read只有落到這些feature上才會(huì)被統(tǒng)計(jì)到，默認(rèn)是“exon”

##Featurecounts結(jié)果說明
Geneid  Chr     Start   End     Strand  Length  ck1.sort.bam
#第一列 Geneid
#第二列 chr，可能對(duì)應(yīng)多個(gè)位置
#第三列 起始位點(diǎn)，可對(duì)應(yīng)多個(gè)位點(diǎn)
#第四列 終止位點(diǎn)，跟前邊一一對(duì)應(yīng)
#第五列 +-鏈信息
#第六列 基因長度
#第七列 樣本表達(dá)信息（最重要的一列）


##提取信息做成數(shù)據(jù)表達(dá)矩陣
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語言統(tǒng)計(jì)比較方便了

3 差異表達(dá)基因統(tǒng)計(jì)

說明文件請(qǐng)參考RNA-seq(6): reads計(jì)數(shù)，合并矩陣并進(jìn)行注釋
計(jì)數(shù)分為三個(gè)水平： gene-level, transcript-level, exon-usage-level
標(biāo)準(zhǔn)化方法： FPKM RPKM TMM TPM
該文使用DEseq2

setwd("./")