RNA-seq分析（二）Alignment

1. STAR （Spliced Transcripts Alignment to a Reference）

1.1 STAR 构建索引

下载基因组及注释文件网站
Saccharomyces cerevisiae (ID 15) - Genome - NCBI (nih.gov)

#download genome 
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz
#download gff
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.gff.gz
#解压 gunzip

正式构建基因组 index

STAR --runThreadN 8 --runMode genomeGenerate --genomeDir index_dir --genomeFastaFiles GCF_000146045.2_R64_genomic.fna --sjdbGTFfile GCF_000146045.2_R64_genomic.gff 

--runThreadN：线程数;
--runMode genomeGenerate： runMode 运行程序模式，默认是比对;
--genomeDir：构建索引输出文件的目录, index_dir；
--genomeFastaFiles：基因组文件, GCF_000146045.2_R64_genomic.fna；
--sjdbGTFfile：基因组注释文件, GCF_000146045.2_R64_genomic.gff；
--sjdbOverhang：reads长度减1

1.2. STAR 比对

建索引的STAR版本，后续比对也用对应的STAR版本，否则报错。

vim star.sh
#!bin/bash
for i in `seq 4 1 9`
do
  STAR --runThreadN 10 --genomeDir ~/rnaseq/1_rnaseq/yeast_genome/index_dir --readFilesCommand gzip -c --readFilesIn ~/rnaseq/SY14/fastq/clean_data/SRR705970${i}_1_val_1.fq.gz ~/rnaseq/SY14/fastq/clean_data/SRR705970${i}_2_val_2.fq.gz --outFileNamePrefix SRR705970${i}_STARdata --outSAMtype BAM SortedByCoordinate
done
qsub -N STAR  -cwd star.sh

如果序列文件是压缩之后的，需要用--readFilesCommand参数指定文件解压缩的方法。gzip -c，--stdout将解压缩的内容输出到标准输出，原文件保持不变；
--readFilesIn：输入fastq文件的路径;
--outFileNamePrefix：输出文件的前缀；
--outSAMtype BAM SortedByCoordinate：输出排序后的bam文件；
outSAMtype该参数有两个字段值，第一个值指定文件类型，取值有SAM和BAM两种；第二个值指定是否排序，取值范围包括Unsorted，SortedByCoordinate；
例如：输出的SRR7059704_STARdataLog.final.out （Log.final.out"里记录了许多比对情况的统计信息）
可参考：STAR:转录组数据比对工具简介 - 云+社区 - 腾讯云 (tencent.com)

2. htseq-count 统计reads counts

HTSeq作为read counts的计数软件，承接的是上游比对软件对于clean data给出的比对结果即bam文件（由sam文件sort得到）。
注意事项：1. HTSeq是对有参考基因组的转录组测序数据进行表达量分析的，其输入文件必须有SAM和GTF文件。

2. -a参数：htseq-count默认会根据mapping的质量值对BAM文件进行过滤，默认值为10, 意味着只有mapping quality > 10的reads才会用来计数。
3. 输入的GTF文件中不能包含可变剪接信息，否则HTSeq会认为每个可变剪接都是单独的基因，导致能比对到多个可变剪接转录本上的reads的计算结果是ambiguous，从而不能计算到基因的count中。
   讲解可参考 【转录组】RNA-seq分析htseq-count的使用 - 喻宇烨 - 博客园 (cnblogs.com)

准备工作

gffread  GCF_000146045.2_R64_genomic.gff -T -o 
#gffread把gff文件转化为gtf, GCF_000146045.2_R64_genomi.gtf
vim samtools_index.sh
#!bin/bash
for i in `seq 4 1 9`
do
  samtools index SRR705970${i}_STARdataAligned.sortedByCoord.out.bam
done
qsub -N samtools_index  -cwd samtools_index.sh
#给bam文件建索引，生成bai文件

正式开始啦~

vim  htseq-count.sh
#!bin/bash
for i in `seq 4 1 9`
do
  htseq-count -f bam -r pos SRR705970${i}_STARdataAligned.sortedByCoord.out.bam ~/rnaseq/1_rnaseq/yeast_genome/GCF_000146045.2_R64_genomi.gtf > SRR705970${i}.count 2 > SRR705970${i}.log
done
qsub -N htseq-count -q g5.q -cwd htseq-count.sh

-f 设置输入文件格式，bam、sam；默认为sam文件；
-r 设置排序方式，name、pos；pos 排序；
例如：">" 标准输出到SRR7059704.count 文件；"2>" 标准误输出到SRR7059704.log文件