參考教程：
數(shù)據(jù)：線蟲轉(zhuǎn)錄組測序PF data（雙端測序）;
流程：
- cutadapt去接頭;
- hisat2建立索引（indexing）;
- hisat2+samtools序列比對（alignment）;
通過index進行比對可加快比對速度。
- FastQC對當前測序數(shù)據(jù)的質(zhì)量進行評估;
- Stringtie對數(shù)據(jù)進行拼接赚导、定量（expression）;
一個基因可能由幾段不相鄰的序列組成锻霎。
- prepDE.py腳本提取基因表達矩陣;
- DESeq2差異分析;
- Metascape路徑分析差異基因富集（Pathway and Process Enrichment Analysis）。
cutadapt去接頭（adapter trimming）
已知接頭序列read1_adapt长已、read2_adapt，需要將read1_adapt 的反向互補序列read1_adapt_REV（在線轉(zhuǎn)換）
cutadapt -a read2_adapt -A read1_adapt_REV -m 20 --pair-filter=both -o out_fq1 -p out_fq2 fq1 fq2
#-m：表示去除接頭后如果read長度小于這個m值就不要了
hisat2建立索引（indexing）
-
hisat2官網(wǎng)
官網(wǎng)提供了一些物種的index文件昼牛，如人類术瓮、小鼠 - 基因注釋：Caenorhabditis_elegans.gtf
- 參考基因組：Caenorhabditis_elegans.WBcel235.dna.toplevel.fa
#加入-ss, -exon，需要消耗上百G的內(nèi)存贰健，不加也可胞四，具體區(qū)別不太清楚
hisat2_extract_splice_sites.py Caenorhabditis_elegans.gtf > Caenorhabditis_elegans.ss
hisat2_extract_exons.py Caenorhabditis_elegans.gtf > Caenorhabditis_elegans.exon
hisat2-build -p 20 Caenorhabditis_elegans.WBcel235.dna.toplevel.fa --ss Caenorhabditis_elegans.ss --exon Caenorhabditis_elegans.exon Caenorhabditis_elegans
#-p 線程數(shù)
#最后一項為輸出的index文件名
hisat2+samtools序列比對（alignment）
hisat2 -x Caenorhabditis_elegans(index_name) -p 20 -1 N1_R1.fastq.gz -2 N1_R2.fastq.gz -S N1.sam
輸出：
26484487 reads; of these:
26484487 (100.00%) were paired; of these:
2355777 (8.89%) aligned concordantly 0 times
23433490 (88.48%) aligned concordantly exactly 1 time
695220 (2.63%) aligned concordantly >1 times
----
2355777 pairs aligned concordantly 0 times; of these:
605991 (25.72%) aligned discordantly 1 time
----
1749786 pairs aligned 0 times concordantly or discordantly; of these:
3499572 mates make up the pairs; of these:
2745181 (78.44%) aligned 0 times
713028 (20.37%) aligned exactly 1 time
41363 (1.18%) aligned >1 times
94.82% overall alignment rate
輸出sam文件，將其轉(zhuǎn)為bam文件（sam的二進制格式）伶椿，并對bam文件進行排序辜伟。
samtools sort -@ 8 -o *.bam *.sam
FastQC
對測序數(shù)據(jù)（去接頭）的質(zhì)量進行評估。
fastqc N1.bam N2.bam N3.bam T1.bam T2.bam T3.bam
Stringtie對數(shù)據(jù)進行拼接脊另、定量（expression）
stringtie -p 8 -G /……/Caenorhabditis_elegans.gtf -e -B -o N1/transcripts.gtf -A N1/gene_abundances.tsv /……/N1.bam
設置參數(shù)-B导狡，輸出的結(jié)果為ballgown所需要的格式，需要python腳本prepDE.py偎痛，才能得到基因表達矩陣旱捧。
使用prepDE.py腳本提取基因表達矩陣
進入ballgown文件夾（由stringtie生成，包含所有樣本）踩麦，將prepDE.py腳本拷貝至當前文件夾廊佩，
cp /……/prepDE.py ./
使用python命令直接運行腳本，注意需要使用python2
python prepDE.py
運行結(jié)果中"gene_count_matrix.csv"是基因表達矩陣（DESeq2的輸入文件）靖榕。
刪除gene_count_matrix.csv中全是0的行（R語言）：
gene_count <- read.csv("gene_count_matrix.csv",stringsAsFactors = F,header = T)
rownames(gene_count) <- gene_count[,1]
gene_count <-gene_count[,-1]
gene_count<-gene_count[which(rowSums(gene_count) > 0),]
DESeq2差異分析
基于負二項式分布标锄，包含標準化（基于負二項式分布）。
參考教程：
# 設置工作目錄
setwd("C:/Users/DELL/Desktop/gene_count/")
# 讀入基因表達值茁计，設定行名為gene_id
gene_count <- read.csv("gene_count_matrix.csv",stringsAsFactors = F)
# 對gene_id一列進行拆分料皇，去除重復名稱
library(stringr)
# 設置空的"gene_count_1"向量，行數(shù)與上面的測序結(jié)果一致
gene_count_1<-rep(NA,nrow(gene_count))
# 利用for循環(huán)星压，對gene_count數(shù)據(jù)框中的重復列進行拆分提取
#將gene_id中|兩側(cè)拆分践剂、提取
for (i in 1:nrow(gene_count)){
gene_count_1[i] <- unlist(str_split(gene_count[i,1], pattern = "\\|"))[1]
}
# 對原數(shù)據(jù)框中的特定序列重新賦值
gene_count$gene_id <- gene_count_1
# 將第一列作為文件的行名
rownames(gene_count) <- gene_count[,1]
gene_count <-gene_count[,-1]
#保存文件至對應目錄
write.csv(gene_count, file = "……/gene_count/gene_count.csv", row.names = TRUE)
#加載DESeq2包
library(DESeq2)
#DESeq2需要三種矩陣，分別為countData表達矩陣,colData樣品信息矩陣及design差異表達矩陣
#countData為表達矩陣即gene_count
#colData為樣品信息矩陣即coldata
#design為差異表達矩陣即批次和條件（對照娜膘、處理）等
#設置condition樣品組別逊脯、重復數(shù)
condition<- factor(c(rep("N", 3), rep("T", 3)), levels = c("N","T"))
condition
#設置樣品信息矩陣colData
colData <- data.frame(row.names = colnames(gene_count), condition)
colData
#在R里面用于構(gòu)建公式對象，~左邊為因變量竣贪，右邊為自變量军洼。
#標準流程：dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition)
#countData為表達矩陣即countdata
#colData為樣品信息矩陣即coldata
#design為差異表達矩陣即批次和條件（對照巩螃、處理）等
#對dds進行標準流程構(gòu)建
dds <- DESeqDataSetFromMatrix(gene_count, colData, design = ~condition)
#過濾低豐度數(shù)據(jù)
dds <- dds[rowSums(counts(dds)) > 1, ]
#或者在構(gòu)建dds之前加上gene_count <- gene_count[apply(gene_count, 1, sum) > 1 , ]
#對原始dds進行normalize
dds <- DESeq(dds)
#顯示dds信息
dds
# 對差異分析結(jié)果進行保存 -------------------------------------------------------------
#使用DESeq2包中的results()函數(shù)，提取差異分析的結(jié)果
#Usage:results(object, contrast, name, .....）
#將提取的差異分析結(jié)果定義為變量"res"
#contrast: 定義誰和誰比較匕争，處理組在前避乏，對照組在后
#提取分析結(jié)果并保存為res
res = results(dds, contrast=c("condition","T","N"))
#對結(jié)果res利用order()函數(shù)按pvalue值進行排序
#order()函數(shù)先對數(shù)值排序，然后返回排序后各數(shù)值的索引甘桑，常用用法：V[order(V)]或者df[order(df$variable),]
#對res進行排序
res = res[order(res$pvalue),]
#顯示res結(jié)果前6行信息
head(res)
#對res矩陣進行總結(jié)拍皮，利用summary命令統(tǒng)計顯示一共多少個genes上調(diào)和下調(diào)
summary(res)
#將差異分析的所有結(jié)果進行輸出保存
write.csv(res, file="……/all_different_genes.csv")
#利用table函數(shù)統(tǒng)計顯著差異基因的數(shù)目
table(res$padj<0.05&abs(res$log2FoldChange>1))
#對具有顯著性差異的結(jié)果進行過濾、提取
#獲取padj小于0.05跑杭，表達倍數(shù)取以2為對數(shù)后的絕對值大于1
#使用subset()函數(shù)過濾需要的結(jié)果至新的變量significant_different_genes_group中
#Usage:subset(x, ...)铆帽，其中x為objects，...為篩選參數(shù)或條件
#對group中數(shù)據(jù)進行過濾德谅、提取
significant_different_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
#使用dim函數(shù)查看該結(jié)果的維度爹橱、規(guī)模
dim(significant_different_genes)
#顯示結(jié)果的前6行信息
head(significant_different_genes)
#對顯著差異基因進行輸出保存
write.csv(significant_different_genes, file = "……/significant_different_genes.csv")
Metascape路徑分析差異基因富集（Pathway and Process Enrichment Analysis）
-
Metascape
step1：輸入差異基因；
上調(diào)基因（up）：log2FoldChange>0
step2：選擇物種女阀；
step3：express analysis（快速分析），custom analysis（自定義分析）屑迂。
