测序数据回来了该怎么办？

相信大家在研究生涯或多或少都会接触到生物信息，以为这是一块很神秘很高深的领域，其实并不难，只要你去看去学去实践，一切都有可能。 本篇主要告诉大家，如果手里有转录组测序的raw data，该怎么做上游分析，下游当然是可以交给我们的R软件去做啦。 1.数据准备 1.1测序数据（reads） 已有fastq文件，ILLUMINA公司的，具体可以查看你手头的测序报告，或一开始的实验设计。 1.2目标物种基因组数据【基因组fa (genome.fa)和注释文件 (gtf/gff3)】 这一步可以从ENSEMBL下载。资源汇总地址：http://ftp.ensembl.org/pub/可以选择最新的 104版本（截止2021/11/1）。

在104版本中选gtf来下载最新的gtf注释文件（步骤同下），选择fasta来选择最新的基因组文件。基因组文件依次选择release-104 > fasta > homo_sapiens > dna 然后在诸多文件中选择*GRCH38.dna.primary_assembly* 服务器下载基因组命令：

wget https://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz --no-check-certificate

复制代码

下载gtf：

wget https://ftp.ensembl.org/pub/current_gtf/homo_sapiens/Homo_sapiens.GRCh38.104.gtf.gz --no-check-certificate

复制代码 2.测序reads分析过程 2.1质控 FastQC 察看数据质量，Trimmomatic 来进行接头，低质量测序reads的过滤与修剪；

#质控
#测序目录下运行，构建fastqc nohup后台命令，-o 输出文件目录
ls /RawData/*gz |xargs -I [] echo 'nohup fastqc [] -o /fastqc &'>fastqc.sh
./fastqc.sh
#将sh和nohup日志搬出至新的文件夹，使文件夹内只剩余zip和html，运行multiqc总结

multiqc ./

#结果得html去浏览器即可总览数据的质量

复制代码

trimmomatic的安装这里不做介绍，目录如命令中/tools/trimmomatic/Trimmomatic-0.39 ps：原始数据命名：T1_1.fq.gz和T1_2.fq.gz,是一个样本双端测序。 trim.sh，在trim目录下运行，这样输出文件方便归纳。

#!/bin/bash



for i in T1 T2 T3 N1 N2 N3

do

java -jar /tools/trimmomatic/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 20 -phred33 /RawData/{i}_1.fq.gz /RawData/{i}_2.fq.gz {i}_1.clean.fastq {i}_1.unpaired.fastq {i}_2.clean.fastq {i}_2.unpaired.fastq ILLUMINACLIP:/tools/trimmomatic/Trimmomatic-0.39/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50
done

复制代码

下面对该软件的参数 进行介绍 PE：双端； threads：线程数； phred编码方式，"新时代"基本都是33而不是64； 当然我们可以简单靠一行代码判断，解压其中一个单端文件N1试试；

head N1.fq | awk '{if(NR%4==0) printf("%s",i; if(i<min) min=

复制代码

返回结果：Phred+33，那么命令中记为33； 接头和引物序列在 TruSeq3-SE 中，第一步 seed 搜索允许2个碱基错配，palindrome 比对分值阈值 30，simple clip 比对分值阈值 10； palindrome 模式允许切除的最短接头序列为 8bp（默认值），palindrome 模式去除与 R1 完全反向互补的 R2（默认去除），上述命令中已省略； ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>:<minAdapterLength>：<keepBothReads>参数中还剩余第五个没有说； 第五个参数keepBothReads：只对 PE 测序的 palindrome clip 模式有效，这个参数很重要，在上图中 D 模式下， R1 和 R2 在去除了接头序列之后剩余的部分是完全反向互补的，默认参数 false，这也就意味着整条去除与 R1 完全反向互补的 R2，当做重复去除掉，但在有些情况下，例如需要用到 paired reads 的 bowtie2 流程，就要将这个参数改为 true，否则会损失一部分 paired reads。这一个上述命令也省略； SLIDINGWINDOW:4:15 #设置 4bp 窗口，碱基平均质量值阈值 15； LEADING是从 reads 的起始端开始切除质量值低于设定的阈值的碱基，直到有一个碱基其质量值达到阈值。 TRAILING是从 reads 的末端开始切除质量值低于设定阈值的碱基，直到有一个碱基质量值达到阈值。 Illumina 平台有些低质量的碱基质量值会被标记为 2 ，所以设置为 3 可以过滤掉这部分低质量碱基。 官方推荐使用 Sliding Window 或 MaxInfo 来代替 LEADING 和 TAILING。（这里未做替换） MINLEN设定一个最短 read 长度，当 reads 经过前面的过滤之后，如果保留下来的长度低于这个阈值，整条 reads 被丢弃。被丢弃的 reads 数会被统计在 Trimmomatic 日志的 dropped reads 中。（也有设置51的) 2.2比对 -x是前面索引解压后的目录，-p线程数，sam目录下运行，-dta是官方推荐后续若使用stringtie组装转录本需要加上的参数

#!/bin/bash



for i in T1 T2 T3 N1 N2 N3

do

hisat2 --dta -x /index/homo_GRCh38.dna.primary -1 /trim/{i}_1.clean.fastq -2 /trim/{i}_2.clean.fastq -S 
done

复制代码

2.3格式转换及排序 -@线程数，bam目录下运行

#!/bin/bash



for i in T1 T2 T3 N1 N2 N3

do

samtools view -bS /sam/{i}.sam | samtools sort -@ 16 -o {i}.sorted.bam
done

复制代码

2.4序列组装 用StringTie对每个样本进行转录本组装 -G前面下载的gtf注释所在目录，stringtiedata目录下运行

#!/bin/bash



for i in T1 T2 T3 N1 N2 N3

do

stringtie /bam/${i}.sorted.bam \

-o /stringtiedata/${i}.gtf \

-p 16 -G /gtf/Homo_sapiens.GRCh38.104.gtf

done

复制代码

获取所有*.gtf 文件名的列表， 并且每个文件名占据一行

ls /stringtiedata/*gtf > Sample_gtf.txt

复制代码

merge目录下运行下述命令，合并这一批样本的gtf

#! /bin/bash



stringtie --merge -p 16 -G /gtf/Homo_sapiens.GRCh38.104.gtf -o merged.gtf Sample_gtf.txt

复制代码

最后在merge目录下，重新用生成的merged.gtf组装每个样本的gtf

#!/bin/bash



for i in T1 T2 T3 N1 N2 N3

do

stringtie -e -B -p 16 -G merged.gtf -o {i}.gtf /bam/{i}.sorted.bam
done

复制代码

2.5count data 提取 准备上述gtf结果文件sample文件 (sample_lst.txt)，格式如下：

Sample1 /merge/T1.gtf

Sample2 /merge/T2.gtf

Sample3 /merge/T3.gtf

Sample4 /merge/N1.gtf

Sample5 /merge/N2.gtf

Sample6 /merge/N3.gtf

复制代码

提取各样品count data

prepDE.py -i sample_lst.txt

复制代码

py文件后续会分享，测序数据大家可以去公共数据库或者尝试自己的数据。 3.差异表达分析 主要就是准备表型文件和上述的基因或转录本count 文件； 表型数据格式如下 (phenodata.csv)：

sample        group

Sample1 Tumor

Sample2 Tumor

Sample3 Tumor

Sample4 Normal

Sample5 Normal

Sample6 Normal

复制代码

采取DESeq2包进行差异分析

library("DESeq2")

countData <- as.matrix(read.csv("gene_count_matrix.csv", row.names="gene_id"))

colData <- read.csv("phenodata.csv", sep="\t", row.names=1)

all(rownames(colData) %in% colnames(countData))

countData <- countData[, rownames(colData)]

all(rownames(colData) == colnames(countData))

dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ group)

dds <- DESeq(dds)

res <- results(dds)

(resOrdered <- res[order(res$padj), ])

复制代码

后续就是自定义fc和padj定义差异基因，以及可视化(火山图、热图之类的)； R中的部分就介绍上面这一步，算是上下游基因差异表达的衔接部分。

代码中需要用到的输入数据：py文件。