我们进行完RNA-seq后，真实返回的测序数据就是这个样子的
https://img.haomeiwen.com/i12544845/e3b832719ef3fe03.png?imageMogr2/auto-orient/strip|imageView2/2/w/1240
我们直接分从其中的rawdata开始分析 因为此次是双端测序，所以每个样本有两个fastq文件。
RNA-seq上游测序的步骤可以概括为 ：
(所有软件均在python3环境下安装)
- 准备参考基因组文件下载及GTF文件下载
- 质控 （fastqc, multiqc）
- 数据过滤 (trim_galore)
- 序列比对 (star)
- 表达定量 (featureCounts)
1. 准备参考基因组文件下载及GTF文件下载:
NCBI下载GTF及参考基因组：
NCBI网址：https://www.ncbi.nlm.nih.gov/genome/?term=mouse
选择getdata
[图片上传失败...(image-e9df51-1610725710244)]
选择download assembly组装好的，
[图片上传失败...(image-76b397-1610725710244)]
参考基因组：
选择REFSEQ
[图片上传失败...(image-e8aa3-1610725710244)]
[图片上传失败...(image-a12dd8-1610725710244)]
这里有GTF文件（gtf.gz 31M）和参考基因组序列(fna.gz 796M)
[图片上传失败...(image-3d3780-1610725710244)]
2. 质控 （fastqc, multiqc）
质量控制
fastqc sample1_R1.fastq.gz #单一样本质控
fastqc sample*gz #批量质控开头为sample，gz结尾的文件
for i in 'ls *gz';do fastqc $i;done #批量质控，遍历该文件夹下所有gz结尾的文件，do fastqc，然后done结束
multiqc时出现了问题：（应该是版本不对）
[图片上传失败...(image-ac0aef-1610725710244)]
解决方法：(multiqc需要安装在python3.7环境下)
[图片上传失败...(image-d75151-1610725710244)]
https://blog.csdn.net/R_python/article/details/105899702
[图片上传失败...(image-b4d663-1610725710244)]
3. 数据过滤 (trim-galore)
注意接头类型（我们用的illumina universal类型）
[图片上传失败...(image-9f6a6d-1610725710244)]
（我们用的illumina universal类型）
[图片上传失败...(image-eff71b-1610725710244)]
实战中使用批量trim-glarore
for i in `seq 91 98`; do trim_galore --paired --q 30 --gzip --fastqc /data/zijun/1.R_project/1.Study_project/3.20201223_ADmouse_RNASeq/rawdata/C${i}_1.fq.gz /data/zijun/1.R_project/1.Study_project/3.20201223_ADmouse_RNASeq/rawdata/C${i}_2.fq.gz -o /data/zijun/1.R_project/1.Study_project/3.20201223_ADmouse_RNASeq/rawdata/trim_out;done
4. 序列比对（star）
实战选择star软件做比对
（1）首先需要对基因组建立索引，建立索引对应的runMode为genomeGenerate, 基本用法如下
# buid index
STAR --runMode genomeGenerate \\
--runThreadN 20 \\
--genomeFastaFiles NCBI_GRCm39_genomic.fna \\
--genomeDir GRCm39_STAR_db \\
--sjdbGTFfile NCBI_GRCm39_genomic.gtf \\
--sjdbOverhang 149
建立索引需要基因组的fasta和gtf文件，通过genomeFastaFiles和sjdbGTFfile这两个参数分别指定；STAR构建索引需要指定一个输出目录，这个目录必须事先创建好，在该目录下，会生成许多文件，所以必须有写权限；runThreadN指定线程数；sjdbOverhang的值默认为100， 在实际设置时，最佳取值为max(read_length) - 1。我们是PE150，所以是149。
在构建索引时，还支持加入intron的区间信息，通过sjdbFileChrStartEnd指定对应的文件，多个文件用逗号分隔，这种格式的文件是由STAR比对产生的，通常用于2-pass比对模式。
官方推荐基因组的fasta采用primary_assembly版本, 不应该包含alt_scaffold和patches。对于human而言，NCBI的链接如下
（2）运行比对
STAR支持fasta/fastq格式的输入文件，如果序列文件是压缩之后的，需要用readFilesCommand参数指定文件解压缩的方法，对于gzip压缩的文件而言，有以下两种下写法
--readFilesCommand zcat
--readFilesCommand gzip -c
比对完成后，会输出许多文件，包含4个类别
- log文件
- sam文件
- bam文件
- 剪切位点文件
每个文件都有事先定义好的名字，当多个样本同时运行时，为了加以区分，可以通过outFileNamePrefix指定输出文件的前缀。前3种类型的文件都比较容易理解，剪切位点文件实际上是根据mapping情况，估算出来的intron区间的信息，默认的文件名称为SJ.out.tab。
默认输出的比对文件为SAM格式，为了节省磁盘空间，方便下游分析，可以通过outSAMtype参数指定输出bam文件，该参数有两个字段值，第一个值指定文件类型， 取值有SAM和BAM两种，第二个值指定是否排序，取值范围包括Unsorted, SortedByCoordinate, 写法如下
--outSAMtype BAM SortedByCoordinate
上述写法输出排序之后的bam文件。
双端数据比对的基本用法如下 (调用了50个线程，reads序列和文件夹，gtf文件必须在同一个文件夹中)
实战中批量比对：
for i in `seq 91 98`
do
STAR \\
--runThreadN 20 \\
--genomeDir GRCm39_STAR_db \\
--readFilesIn C${i}_1_val_1.fq.gz C${i}_2_val_2.fq.gz \\
--readFilesCommand zcat \\
--sjdbGTFfile NCBI_GRCm39_genomic.gtf \\
--sjdbOverhang 149 \\
--outFileNamePrefix 03_align_out/sample_C${i}_ \\
--outBAMsortingThreadN 5 \\
--outSAMtype BAM SortedByCoordinate \\
--quantMode TranscriptomeSAM GeneCounts
done
[图片上传失败...(image-f27a0c-1610725710243)]
出现错误，添加命令 --outBAMsortingThreadN 5， 参考
运行STAR时遇到的错误https://www.jianshu.com/p/e2a149d9e615
[图片上传失败...(image-dcf56f-1610725710243)]
6. 数据定量（featureCounts）
单一文件测试：
featureCounts -p -a /data/zijun/1.R_project/1.Study_project/rawdata/ref/NCBI_GRCm39_genomic.gtf -o \\
out_counts.txt -T 50 -t exon -g gene_id sample_C91_Aligned.sortedByCoord.out.bam
批量测试：
for i in `seq 91 98`
do
featureCounts -p -a /data/zijun/1.R_project/1.Study_project/rawdata/ref/NCBI_GRCm39_genomic.gtf -o \\
out_counts.txt -T 50 -t exon -g gene_id sample_H${i}_Aligned.sortedByCoord.out.bam
done
# 使用multiqc对featureCounts统计结果进行可视化。
multiqc out_counts.txt.summary
# 只保留all.id.txt文件的【基因名】和【样本counts】
cat all.id.txt | cut -f1,7- > counts.txt
[图片上传失败...(image-1784bd-1610725710243)]
[图片上传失败...(image-21f12b-1610725710243)]
所有样本：
featureCounts -p -a /data/zijun/1.R_project/1.Study_project/rawdata/ref/NCBI_GRCm39_genomic.gtf -o \
out_ALL_counts.txt -T 50 -t exon -g gene_id sample_*_Aligned.sortedByCoord.out.bam
网友评论