1.构建参考基因组的索引
1.1以小鼠的基因组为例,首先下载参考基因组与基因组的注释文件
wget http://ftp.ensembl.org/pub/release-102/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
wget http://ftp.ensembl.org/pub/release-102/gtf/mus_musculus/Mus_musculus.GRCm38.102.gtf.gz
#不能使用压缩格式的基因组和注释,先解压
gunzip Mus_musculus.GRCm38.102.gtf.gz Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
1.2调用STAR构建索引
STAR --runThreadN 16 --runMode genomeGenerate --genomeDir /apps/users/zebrafish/PRJNA285844/reference_genome --genomeFastaFiles Mus_musculus.GRCm38.dna.primary_assembly.fa --sjdbGTFfile Mus_musculus.GRCm38.102.gtf --sjdbOverhang 49
--runThreadN:使用的线程数
--runMode:这里填写‘genomeGenerate’,意思是构建基因组索引
--genomeDir:产生索引的存放目录
--genomeFastaFiles:参考基因组的文件,另外可以指定一个或者多个包含基因组的.fa文件,没有这方面的应用,所以不清楚什么实际问题可以用这个选项来解决。一般来说模式生物只用一个参考基因组就行,不用特意制定每个染色体。
--sjdbOverhang :填写测序长度-1,比如这组数据的为SE50测序所以填写49
1.3mapping reads to genome
STAR --runThreadN 16 --genomeDir /apps/users/zebrafish/PRJNA285844/reference_genome --readFilesCommand zcat --readFilesIn clean_SRR2050895.fastq.gz --twopassMode Basic --outSAMstrandField intronMotif --outSAMtype BAM SortedByCoordinate --outFileNamePrefix /apps/users/zebrafish/PRJNA285844/simply_filter/align/SRR2050895
#比对过程中可以修改gtf的路径和sjdbOverhang 的大小
STAR比对链特异性库不需要别的只需要加个标签--outSAMstrandField intronMotif
--twopassMode Basic 比对两遍更准确。
--outSAMtype BAM SortedByCoordinate 输出文件为sort后的bam
1.4 基于count统计表达量
使用featurecount这个软件
featureCounts -T 16 -s 1 -M -t exon -g gene_id --extraAttributes gene_name -a /apps/users/zebrafish/PRJNA285844/reference_genome/Mus_musculus.GRCm38.102.gtf -o SRR2050895 SRR2050895Aligned.sortedByCoord.out.bam
-T:16线程
-s:1 这组数据是正义链的链特异性库
-M:对于多比对的reads都进行统计
-t : 想要统计的feature 的名称, 取值范围是gtf 文件中的第3列的值,默认是exon
-g : 想要统计的meta-feature 的名称,取值范围参考gtf 第9列注释信息,gtf 的第9列为 key=value 的格式, -g 参数可能的取值就是所有的key, 默认值是gene_id
-a:参考的基因组注释文件
-o:输出文件的prefix
然后截取部分需要的col
cut -f 1,6,7,8 SRR2050895 |grep -v '^#' >SRR2050895_Counts.txt