1.索引建立
hisat2-build 基因组.fasta 输出结果文件位置
备选：根据比对的效果，可以另选多种比对方式，也可magicblast，那此时建立索引及比对，见
https://www.jianshu.com/p/ece803eb552b
2.与基因组进行比对
nohup hisat2 -x 基因组索引的前缀 -1 xx.fq1.gz -2 xx.fq2.gz -S 输出文件/xx.sam > 输出文件/xx.log 2>&1 &
Hisat2输出结果
比对之后会输出如下结果，解读一下就是全部数据都是100%的，30.36%的配对数据一次都没有比对，26.31%的数据比是唯一比对，43.33%是多个比对。然后如果不按照顺序来，有1.22%的比对。之后把剩下的部分用单端数据进行比对的话，65.30%数据没比对上，15.06%的数据比对一次，19.64%比对超过一次。总体来看有80.42%的比对。
60021748 reads; of these:
60021748 (100.00%) were paired; of these:
18220061 (30.36%) aligned concordantly 0 times
15792611 (26.31%) aligned concordantly exactly 1 time
26009076 (43.33%) aligned concordantly >1 times
----
18220061 pairs aligned concordantly 0 times; of these:
223096 (1.22%) aligned discordantly 1 time
----
17996965 pairs aligned 0 times concordantly or discordantly; of these:
35993930 mates make up the pairs; of these:
23504789 (65.30%) aligned 0 times
5419179 (15.06%) aligned exactly 1 time
7069962 (19.64%) aligned >1 times
80.42% overall alignment rate
2. SAMtools三板斧
3.samtools 将 sam 转成bam
samtools view -S xx.sam -b > xx.bam
4.samtools将得到的bam文件分类排序
samtools sort xx.bam -o xx.sorted.bam &
5.featureCounts 对CDS进行计数
准备基因组注释文件.gtf
featureCounts --primary -M -T 64 -t CDS -g gene_id -a $gtf -o featurecount.txt *.sort.bam
可以不看此，这是针对每一个组分别进行feature定量：
for id in $(ls ./*sorted.bam | sed 's/.sorted.bam//');do echo "featureCounts --primary -M -T 10 -t CDS -g transcript_id -a /nfs_genome/home/zhaoshanzhong/zs/data/GIPL_genome/OMTSEQ20210416.4_BAPL.genome_pan.result_20211105/Pangenome_ALLmerge.fa.gtf -o $id.featurecount.txt $id.sorted.bam";done > feature.sh
6.DEGseq.R 脚本
vim 生成design 与 comparsion文件（design文件中必须将featurecounts.txt表格中所有的分组信息列备完全，comparison文件可以针对自己所要比较的两组单独生成，不然导致生成的差异基因数目不准）
样品名生成：ll *sorted.bam | awk '{print $9}' | awk -F '.' '{print $1 "\t" $0}'
然后DEGseq.R -c featurecount.txt -C xxcomprsion -o xxvsxx -d design
7.富集分析
①diamond blastp --id 40 --query-cover 65 -e 1e-8 -k 1 -d /nfs_genome/BioDB/eukaryotes.dmnd -q ../output/annot.faa -o anno.m8
diamond blastp --id 40 --query-cover 65 -e 1e-8 -k 30 -d /nfs_genome/BioDB/eukaryotes.dmnd -q ../output/annot.faa -o anno.m8
②blast2kegg.pl folder
利用.m8.kegg文件生成gene2ko文件
awk '$9~/^K/ {print $1" "$9}' anno_lf.m8.kegg | sed "s/\(K[0-9]\+\),.*/\1/" >lifo_gene2ko.tsv
生成的文件中，xxgene2map.csv 以及xxgene2ko.tsv 需要进行对第一列的蛋白数据进行替换修改掉后面蛋白版本号，常用vim中 ：%s/xx//g 进行替换修改，此时vim的常见替换修改方法见：
https://blog.csdn.net/weixin_30471013/article/details/116973791
https://blog.csdn.net/legend050709/article/details/120856828
https://blog.csdn.net/windyf2013/article/details/107229812?
修改好gene2map 和gene2ko表格后即可进行下述命令：
KEGGenrich.r
KEGGenrich.R -k anno_lf2.m8_gene2map.csv -g lifo2_gene2ko.tsv -s ..DEGseq.R生成的tsv文件 -o kegg -c featurecounts生成的文件
③blast2go
该分析需要付费软件，单独联系老板，另外需要自己先利用基因组的蛋白信息生成xml格式的nr注释的比对文件，方法如下：
diamond blastp --sensitive -e 1e-4 -k 30 -d /nfs_genome/BioDB/nr.202108.dmnd -q Pangenome_ALLmerge.faa -o xx.nr.anno --outfmt 5
下一步需要blast2go 界面版本在线注释，获得了注释的文件
blast2go_go_propagation.txt后进行下一步
GOenrich.r
GOenrich.r -g xxblast2go_go_propagation.txt -s degxx.tsv（差异基因表格） -o degxxGO
