转录组分析（最终）

cpu=4
###数据质控###
##质控查看，根据质控结果决定剪去前多少bp##
fastqc -o ./ -t $cpu 随便看一个.fq.gz
##过滤低质量数据。一般同一批下机数据都差不多，可以写个for什么的全用同样的参数过滤##
fastp -q 30 -w $cpu -f 15 -F 15 -i R1.fq.gz -I R2.fq.gz -o R1.clean.fq.gz -O .R2.clean.fq.gz 
#-q 质量设定，-q 30 即设置过滤低于Q30的数据。不过fastp的过滤不怎么狠，还是会有少量低于设定质量值的被保留；
#-f 剪去R1前多少bp，-F 剪R2；
#-i 输入R1，-I 输入R2；
#-o 输出R1，-O 输出R2；
###转录组数据比对与组装###
datapath=yourpath
genomefa=yourgenome.fa
genomegff=yourgenome.gff3
gffread $genomegff -T -o genome.gtf
#获取剪接位点、外显子信息（可做可不做），建立基因组索引#
extract_splice_sites.py genome.gtf > genome.ss
extract_exons.py genome.gtf > genome.exon
hisat2-build --ss genome.ss --exon genome.exon $genomefa genome
##与参考基因组比对，组装转录本（包括计算read counts、FPKM、TPM）。可以一套跑下来，但需要关注mapping rate，起码要有7、80%，不能太低##
for i in `ls *.R1.clean.fq.gz`
do
 i=${i/.R1.clean.fq.gz/}
 hisat2 -p $cpu --dta -x genome -1 $i'.R1.clean.fq.gz' -2 $i'.R2.clean.fq.gz' -S $i'.sam'
 samtools sort -@ $cpu -o $i'.bam' $i'.sam'
 rm $i'.sam'
 #计算相对表达量FPKM和TPM#
 stringtie -e -p $cpu -G genome.gtf -A $i'.tab' $i'.bam'
 #生成gtf文件用于后续计算read counts#
 stringtie -e -p $cpu -G genome.gtf -o $i'.stringtie.gtf' $i'.bam'
done