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

生信喵

相信大家在研究生涯或多或少都会接触到生物信息，以为这是一块很神秘很高深的领域，其实并不难，只要你去看去学去实践，一切都有可能。
本篇主要告诉大家，如果手里有转录组测序的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*
      服务器下载基因组命令：

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

复制代码

     下载gtf：

1. 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的过滤与修剪；

1. #质控
   
2. #测序目录下运行，构建fastqc nohup后台命令，-o 输出文件目录
   
3. ls /RawData/*gz |xargs -I [] echo 'nohup fastqc [] -o /fastqc &'>fastqc.sh
   
4. ./fastqc.sh
   
5. #将sh和nohup日志搬出至新的文件夹，使文件夹内只剩余zip和html，运行multiqc总结
   
6. multiqc ./
   
7. #结果得html去浏览器即可总览数据的质量

复制代码

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

1. #!/bin/bash
   
2. 
   
3. for i in T1 T2 T3 N1 N2 N3
   
4. do
   
5. 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
   
6. done

复制代码

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

1. 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组装转录本需要加上的参数

1. #!/bin/bash
   
2. 
   
3. for i in T1 T2 T3 N1 N2 N3
   
4. do
   
5. 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
   
6. done

复制代码

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

1. #!/bin/bash
   
2. 
   
3. for i in T1 T2 T3 N1 N2 N3
   
4. do
   
5. samtools view -bS /sam/${i}.sam | samtools sort -@ 16 -o ${i}.sorted.bam
   
6. done

复制代码

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

1. #!/bin/bash
   
2. 
   
3. for i in T1 T2 T3 N1 N2 N3
   
4. do
   
5. stringtie /bam/${i}.sorted.bam \
   
6. -o /stringtiedata/${i}.gtf \
   
7. -p 16 -G /gtf/Homo_sapiens.GRCh38.104.gtf
   
8. done

复制代码

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

1. ls /stringtiedata/*gtf > Sample_gtf.txt

复制代码

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

1. #! /bin/bash
   
2. 
   
3. stringtie --merge -p 16 -G /gtf/Homo_sapiens.GRCh38.104.gtf -o merged.gtf Sample_gtf.txt
   

复制代码

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

1. #!/bin/bash
   
2. 
   
3. for i in T1 T2 T3 N1 N2 N3
   
4. do
   
5. stringtie -e -B -p 16 -G merged.gtf -o ${i}.gtf /bam/${i}.sorted.bam
   
6. done

复制代码

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

1. Sample1 /merge/T1.gtf
   
2. Sample2 /merge/T2.gtf
   
3. Sample3 /merge/T3.gtf
   
4. Sample4 /merge/N1.gtf
   
5. Sample5 /merge/N2.gtf
   
6. Sample6 /merge/N3.gtf

复制代码

     提取各样品count data

1. prepDE.py -i sample_lst.txt

复制代码

     py文件后续会分享，测序数据大家可以去公共数据库或者尝试自己的数据。


3.差异表达分析
      主要就是准备表型文件和上述的基因或转录本count 文件；
      表型数据格式如下 (phenodata.csv)：

1. sample        group
   
2. Sample1 Tumor
   
3. Sample2 Tumor
   
4. Sample3 Tumor
   
5. Sample4 Normal
   
6. Sample5 Normal
   
7. Sample6 Normal

复制代码

     采取DESeq2包进行差异分析

1. library("DESeq2")
   
2. countData <- as.matrix(read.csv("gene_count_matrix.csv", row.names="gene_id"))
   
3. colData <- read.csv("phenodata.csv", sep="\t", row.names=1)
   
4. all(rownames(colData) %in% colnames(countData))
   
5. countData <- countData[, rownames(colData)]
   
6. all(rownames(colData) == colnames(countData))
   
7. dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ group)
   
8. dds <- DESeq(dds)
   
9. res <- results(dds)
   
10. (resOrdered <- res[order(res$padj), ])

复制代码

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

生信喵

prepDE.py文件呢，不是我自己弄的哈，是官方还是哪位大佬写的。我下面会分享出来

生信喵

代码中需要用到的输入数据：临床信息和转录组数据。获取示例数据请在公众号"生信喵实验柴"后台回复“20211101”。