生信初入门——转录组数据分析（一）

关于生信处理在下啥都不会，老师让我先从RNA-seq入手学习（听说这个最适合新手 0_0）。

手里有三个RNA-seq的双端测序数据：100cell_PBMC、1cell_PBMC以及对照组的scRNA_PBMC，均无重复。找出前两个实验组分别相对于对照组的差异表达基因。

流程大致为：对测序数据进行质控（linux环境，以下亦是）——将质控好的测序数据比对到基因组上——使用featurecount对比对结果进行counts计数——使用edgeR对counts结果进行差异表达分析（R环境）。

一、 质控

质控这块主要使用的软件为fastp（对数据质量进行操作）和fastqc（查看数据质量如何）。

1. linux环境中进入数据目录，进行以下质控操作

for i in *gz; do fastp -i ${i} -o 输出路径/${i} -h 路径/${i}.html -j 输出路径/${i}.json --thread=4; done

主要使用了一个循环（还不太会写，凑合用了 0_0），使用fastp默认参数将测序数据作为单端进行处理，主要进行低质量reads的过滤、接头剪切等操作。

图1

由于测序过程中随机引物的存在，会导致reads的前端出现杂乱，如下图1，接下来对自动化质量处理后的测序reads的前端进行切除，根据图中的杂乱，主要是对前20bp进行切除。

for i in *gz; do fastp -i ${i} -A -f 20 -o cut_${i} -h cut_${i}.html -j cut_${i}.json --thread=4; done

切后的结果如图2：

图2

（这两部应该是能合并一步处理的，但是不知道为什么我直接一步操作的话，上图波动特别大，就分开来做了 0_0）

2. 将质控后的R1与R2合并为一个文件，后续的所有操作作为单端序列reads进行处理（这步是听说是单细胞的信息较少，将测序数据作为双端处理会丢失些信息，我也没搞太清楚0_0）

以100cell为例：

cat PBMC_100cell_clean_R1.fastq.gz PBMC_100cell_clean_R2.fastq.gz > PBMC_100cell_clean.fastq.gz

（这一步后可以使用fastqc软件对处理后的reads数据进行质控结果的查看：

for file in PBMC*gz; do fastqc -t 4 ${file}; done

也可以不用这部操作，因为fastp中已经包含质量结果的展示了。）

二、 比对

1. 首先将质控后的reads使用hisat2（或STAR）比对到基因组上

hisat2：

for file in PBMC*gz; do hisat2 --dta -p 4 -t -x 参考基因组索引文件目录 -U ${file} -S ${file}.sam; done

STAR：

STAR --runThreadN 16 --readFilesCommandzcat --genomeDir 参考基因组索引文件目录 --readFilesIn PBMC_100cell_clean.fastq.gz --outSAMtype BAM Unsorted --outFileNamePrefix 输出文件前缀

2. 接着将比对结果使用samtools软件进行排序并转换为bam（二进制，占内存小，后面的处理也为bam类型输入）文件类型。

for file in cut*sam; do samtools sort -@ 4 -o ${file}.bam ${file}; done

三、 counts计数

使用subread软件中的featureCounts模块将三个样本得到的比对文件（bam后缀）同时进行counts计数：

featureCounts -T 4 -a 基因组注释文件目录/gencode.v41.annotation.gtf.gz -o total_counts.txt *bam

这样就得到一个包含有三个样本counts的txt文件如图3，后续的操作就是使用该counts文件在R环境中进行差异表达的操作。

图3

(菜鸟一枚，大佬指教，万分感谢 ^_^)

待续 0_0 ...