RNA-seq分析

一、获取数据

1.1 直接下载fastq文件，或者"sratools"获取sra文件后转fastq，下载list.txt后调用prefetch

$ cd /mnt/f/file/sra_gse117489/
$ prefetch --option-file SRR_Acc_List.txt

1.2 将sra文件整合后转fastq文件

$ for I in `seq 28 50`; do fastq-dump --gzip –split-3 -O ./fastq/ -A SRR3910${I}.sra; done

1.3 "fastq-dump"参数说明

$ fastq-dump \ 
--gzip \ # 输出文件为压缩包
–split-3 \ # 若为双端测序，输出"fastq_1.gz"、"fastq_2.gz"两个文件
-O ./fastq/ \ # 输出文件位置
-A SRR391033.sra # 输入文件

二、"trim_galore"编辑

$ conda activate rnaseq
$ cd /mnt/f/file/gse117489/gsm3301893_star/trim_galore_result/
$ trim_galore \
--output_dir ./ \
--length 75 \
--quality 20 \
--stringency 5 ../967_968sra_data.fastq

三、"fastqc"质控

$ cd /mnt/f/file/gse117489/gsm3301893_star/trim_galore_result/fastq_result_for_new/ # 定位
$ sudo /mnt/f/biosoft/fastqc/FastQC/fastqc -o ./ ../967_968sra_data_trimmed.fq # 运行fasqc

四、"star"比对基因组

4.1 构建索引

$ conda activate rnaseq # 激活安装有star的conda环境
$ cd /mnt/f/file/genome/ # 定位到基因组文件位
$ mkdir star_index # 创建索引文件夹

$ STAR --runThreadN 8 --runMode genomeGenerate --genomeDir ./star_index --genomeFastaFiles ./GRCh38.p13.genome.fa --sjdbGTFfile ./gencode.v39.chr_patch_hapl_scaff.annotation.gtf --sjdbOverhang 100 --limitGenomeGenerateRAM 13966723168 --genomeSAsparseD 2 --limitSjdbInsertNsj 1000000

4.2 "star"构建索引参数说明

$ STAR \
--runThreadN 8 \ # 线程数
--runMode genomeGenerate \ # 模式
--genomeDir ./star_index \ # 输出索引位置
--genomeFastaFiles ./GRCh38.p13.genome.fa \ # 基因组文件
--sjdbGTFfile ./gencode.v39.chr_patch_hapl_scaff.annotation.gtf \ # 基因组注释文件
--sjdbOverhang 100 \
--limitGenomeGenerateRAM 13966723168 \ # 限制内存，否则内存被占满后将报错
--genomeSAsparseD 2 \
--limitSjdbInsertNsj 1000000 # 限制？，不加会报错

4.3 "star"比对命令

$ STAR --runThreadN 8 --runMode alignReads --quantMode TranscriptomeSAM GeneCounts --outSAMtype BAM Unsorted --outSAMunmapped None --genomeDir /mnt/f/file/genome/star_index/ --readFilesIn ../SRR7550967_sra_data.fasta ../SRR7550968_sra_data.fasta --outFileNamePrefix gem33

4.3 "star"比对命令参数说明

STAR \
--runThreadN 8 \ # 线程数
--runMode alignReads \ # 模式
--quantMode TranscriptomeSAM GeneCounts \ #
--outSAMtype BAM Unsorted \
--outSAMunmapped None \
--genomeDir /mnt/f/file/genome/star_index/ \ # 基因组索引位置
--readFilesIn \ # 输入文件
--outFileNamePrefix gem33 \ # 输出文件前缀

五、"htseq-count"计数

5.1 "htseq"命令

$ conda activte py2
$ htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m intersection-nonempty htseq_ceshi.bam /mnt/f/file/genome/gencode.v39.chr_patch_hapl_scaff.annotation.gtf >htseq_counts.csv

5.2 "htseq-count"参数说明

$ conda activte py2 # 激活安装有htseq的conda环境
$ htseq-count \
-f bam \ # 输入文件格式
-r name \ # 按名字排序
-s no \
-a 10 \ # 滤过低质量
-t exon \
-i gene_id \
-m intersection-nonempty \
/mnt/f/file/gse117489/gsm3301893_star/star_result_ceshi/htseq_ceshi.bam \ # 输入文件
/mnt/f/file/genome/gencode.v39.chr_patch_hapl_scaff.annotation.gtf \ # 基因注释文件（gtf or gff）
>htseq_counts.txt # 输出 或者 -0 输出文件