将 gff 文件转成 gtf (featurecounts需要使用gtf文件)

gffread coreset.gff -T -o amur_ide.gtf
# -o    write the output records into <outfile> instead of stdout
#-T    main output will be GTF instead of GFF3

构建参考基因组的索引文件

hisat2-build -p 8 genome.fa amur_ide

hisat2批量比对

for i in 39 40 41 42 43 44
do 
nohup hisat2 -x /home/genome_index/amur_ide -1 SRR75089${i}_1.fq -2 SRR75089${i}_2.fq | samtools view -S -b > xx.bam &
done

bam文件排序

samtools sort XX.bam -o xxx_sorted.bam

featurecounts 定量

for i in 39 40 41 42 43 44
do 
nohup featureCounts -p -a /home/jiamj/analysis/ref/TAIR10.gtf -o ${i}_counts.txt /home/jiamj/analysis/clean/${i}_sorted.bam &
done

-p If specified, libraries are assumed to contain paired-end reads. For any library that contains paired-end reads, the 'countReadPairs' parameter controls if read pairs or reads should be counted

结果包含有 geneid，染色体位置，基因起始结束的位置以及基因的 count 数

image.png

featureCounts进行fpkm标准化

countdata <- read.csv("countdata.csv")
#countdata.csv是提取了上一步的counts数据以及gene length
rownames(countdata) <- countdata[,1]
countdata <- countdata[,-1]
kb <- countdata$length / 1000
count <- countdata[,1:8]
rpk <- count / kb
tpm <- t(t(rpk)/colSums(rpk) * 1000000)
fpkm <- t(t(rpk)/colSums(count) * 10^6) 
#想计算数据框中每列的总和，使用colSums函数。
write.table(fpkm,file="eight_tissues_fpkm.xls",sep="\t",quote = F)