转录组分析流程：比对（有参）及统计Counts矩阵

1. 质控

fastqc *
multiqc *
trimmomatic_run.sh #去掉前9个碱基

trimmomatic_run.sh

#! /bin/bash
#history:
#	Gossie V1 20190101

#用trimmomatic对进行测序数据质量控制
#输入文件为包含样本名称的单列文件

#线程
t=1
# 模式PE（双端）或者SE（单端）
mode='PE'

#以下参数适情况修改、添加。
#ILLUMINACLIP模式
#接头文件，fasta格式。注意：不同的命名格式有差别。
fastaWithAdapterEtc='/Trimmomatic-0.38/adapters/TruSeq3-PE-2.fa'
#第一步搜索时允许搞错错配碱基个数，通常为1-2.
seed_mismatches=2
#双端测序的palindrome模式下，连上接头的两条reads的比对阈值
palindrome_clip_threshold=30
#切除的接头序列的最低比对分值，通常为7-15.
simple_clip_threshold=10
#仅对palindrome模式有效，双端测序palindrome模式下，可以切除的接头的最短长度。默认是8。但，可以设置为1.
minAdapterlength=8
#仅对palindrome模式有效。默认为false。
keepBothReads=true

# SLIDINGWINDOW:滑窗模式。从5‘端开始滑动，切除平均质量小于阈值的片段。
#窗口大小
windowSize=0
#质量
requiredQuality=0

#LEADING:切除起位置小于阈值的碱基
lead_quality=0
#TRAILING:切除末尾位置小于阈值的碱基
trail_quality=0
#CROP:从末尾切除read部分碱基，使reads达到指定长度.
crop_length=150
#HEADCROP:切除起始位置特定长度的碱基
headcrop=9
#MINLEN:舍弃小于给定长度的
minlen=50
#AVGQUAL：舍弃平均质量低于给定水平的read
avgqual=20

while read id
do
nohup trimmomatic PE -threads $t ${id}_1.fastq.gz ${id}_2.fastq.gz -baseout ${id}.fastq.gz HEADCROP:$headcrop ILLUMINACLIP:$fastaWithAdapterEtc:$seed_mismatches:$palindrome_clip_threshold:$simple_clip_threshold:$minAdapterlength:$keepBothReads MINLEN:$minlen & 
done < $1

2. 比对

比对软件：hisat2

gffread annotation.gff3 -T -o annotation.gtf #hisat2需要GTF格式的注释文件，注意看一下有没有外显子信息(exon)，有的gff3可能不提供，会影响比对reads的统计
hisat2_extract_splice_sites.py .UTX.gtf >UTX.gene.ss #剪接位点文件

hisat2_extract_exons.py UTX.gene.gtf >UTX.gene.exon #外显子位点文件

sh hisat2.sh #samtools version 1.21

hisat2.sh
需要 1.9以上版本的samtools用于排序和转换bam

#! /bin/bash
# 输入文件：sra id
#index文件
index=/hisat2_index/UTX
while read id
do
#--dta 输出比对情况
nohup hisat2 -p 2 --dta -x $index -1 ${id}_1P.fastq.gz -2 ${id}_2P.fastq.gz -S $id.hisat2.sam >$id.hisat2.log 2>&1 && samtools sort -@ 10 -o $id.sorted.bam $id.hisat2.sam && rm $id.hisat2.sam &
done < $1 

3. 计算count矩阵

软件：featureCounts

nohup featureCounts -p -T 64 -t exon -g gene_id -a UTX.gene.gtf -o counts.featurecounts.txt *.sorted.bam >featurecounts.log 2>&1 &
# -p：双端；单端不用加参数
# -T：线程，1~64
# -t：feature type
# -g：meta-feature 类型，若一个reads或fragment比对到多个feature（如，exon),如果这些exon属于一个meta-feature（如，gene)，则只计数一次
# -a：注释文件，最好用GTF；因为GFF3的ID注释包含gene、mRNA、exon、UTR等
# -o：输出文件
# 可对多个bam同时统计