测序数据回来了该怎么办？
本帖最后由 生信喵 于 2021-11-1 10:59 编辑相信大家在研究生涯或多或少都会接触到生物信息，以为这是一块很神秘很高深的领域，其实并不难，只要你去看去学去实践，一切都有可能。
本篇主要告诉大家，如果手里有转录组测序的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",$0);}' |od -A n -t u1 | awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding!";}' 返回结果：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 ${i}.sam -p 16
done2.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
done2.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
done2.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) 后续就是自定义fc和padj定义差异基因，以及可视化(火山图、热图之类的)；
R中的部分就介绍上面这一步，算是上下游基因差异表达的衔接部分。
prepDE.py文件呢，不是我自己弄的哈，是官方还是哪位大佬写的。我下面会分享出来 代码中需要用到的输入数据：临床信息和转录组数据。获取示例数据请在公众号"生信喵实验柴"后台回复“20211101”。
页:
[1]