生物信息-转录组分析流程

如何下载NCBI中的SRA数据

下载我们需要用到一个sra-tools,可以建立新的conda环境用conda安装

conda install sra-tools

#先找到SRA database中的基因（SRA_accessionList.csv）

#批量下载基因
awk '{print "prefetch" $1 "&"}' SRA_accessionList.csv > run_prefetch.sh
#利用awk生成代码并保存再shell文件中

#将sra转换为fastq
fastq-dump xxx.sra

#下载参考基因组
genome.fasta
genome.gtf

#质控
fastqc xxx.fastq
trimmomatic PE  <input1.fastq> <input2.fastq> <output1_paired.fastq> <output1_unpaired.fastq> <output2_paired.fastq> <output2_unpaired.fastq> ILLUMINACLIP:<adapter.fa>:2:30:10 LEADING:15 TRAILING:15 SLIDINGWINDOW:4:15 MINLEN:36

#比对到参考基因组
hisat2-build genome.fasta genome 1>hisat2-build.log2>&1

SNP挖掘

质控

质控是非常重要的一步，通常我们可以使用fastqc + trimmomatic进行分析 + 质量值过滤两个工具来实现，也可以使用我们最新的国产软件 fastp 来进行质控和过滤一步到位。

fastqc xxx.fastq.gz #生成单个报告
multiqc ./#生成组合报告

双末端数据选PE，变量是参考的
trimmomatic PE  <input1.fastq> <input2.fastq> <output1_paired.fastq> <output1_unpaired.fastq> <output2_paired.fastq> <output2_unpaired.fastq> ILLUMINACLIP:<adapter.fa>:2:30:10 LEADING:15 TRAILING:15 SLIDINGWINDOW:4:15 MINLEN:36

以上是一个使用Trimmomatic进行双端测序样本修剪的示例命令，包含了一些常用的参数设置：

请根据具体情况进行以下替换和调整：

<input1.fastq>和<input2.fastq>：替换为您的双端测序样本的输入文件路径。
<output1_paired.fastq>``<output1_unpaired.fastq>：表示经过处理后第一个样本的配对 reads 和未配对 reads 输出的文件
<output2_paired.fastq>``<output2_unpaired.fastq>：表示经过处理后第二个样本的配对 reads 和未配对 reads 输出的文件
<adapter.fa>：这里提供特定的adapter的fasta文件，2 表示最短的重复匹配长度为 2，30 表示最大允许的不匹配的错误率为 30，10 表示在适配器和序列之间允许的最大差异为 10
参数设置如LEADING:15和TRAILING:15表示移除序列开头和结尾低质量碱基的阈值
SLIDINGWINDOW:4:15指定了滑动窗口的大小和最小平均质量要求，表示滑动窗口的大小为4，过滤质量值的阈值为15
MINLEN:36设置修剪后序列的最小长度
我们可以通过fastqc提供的样品整体质量值分布来确定阈值，根据数据量的大小来确定SLIDINGWINDOW的大小和MINLEN

此外，我们可以使用fastp来进行过滤和质控，这将节省我们的精力来手动设置参数

fastp -i input1.fastq -I input2.fastq -o output/filtered1.fastq -O output/filtered2.fastq

其中，input1.fastq 和 input2.fastq 是你的输入双末端数据文件，output/filtered1.fastq 和 output/filtered2.fastq 是过滤后的输出文件路径。

如果我们这样写的话，fastp会自动对 reads 进行质量修剪，去除低质量的碱基，去除含有接头序列（adapter）的 reads，过滤掉长度过短的 reads，过滤掉头尾质量过低的 reads，这是一定程度上的智能控制，当然我们也可以进一步自行设置参数，以达到我们自己的实验目的。

TIPS:可选的参数：你还可以使用一些可选参数来指定更多的质控过滤设置，例如：

-q <quality>：设置质量阈值，低于此阈值的碱基将被过滤掉。
-u <length>：设置过滤掉长度低于指定值的读序。
-x：执行最小长度过滤，将两个 reads 中较短的一个过滤掉。

根据你的需求，可以根据 fastp 文档中的说明，调整这些参数。

请确保在运行命令之前将上述命令中的文件名 input1.fastq、input2.fastq 和输出路径 output/filtered1.fastq、output/filtered2.fastq 替换为实际文件名和路径。

运行完以上命令后，fastp 将会对输入的双末端数据进行质量控制和过滤，并将过滤后的结果保存到指定的输出文件中。

将测序数据比对到参考基因组,并生成比对结果sam文件

最基础的比对算法是bwa,BWA 是一个常用的基因组比对工具，它使用 Burrows-Wheeler Transform（BWT）索引算法来实现快速且高效的比对。使用它通常需要先构建索引，然后比对到参考基因组，参考基因组的获得可以通过组装基因组或者在open resource下载得到。

bwa index ref.fa#构建索引，增加比对速度

bwa mem -t 10 ref.fa xxx_1.fastq.gz xxx_2.fastq.gz > map.sam#比对

在使用 BWA（Burrows-Wheeler Aligner）等比对工具对参考基因组进行比对后会生成以下文件：

在比对过程中，BWA 会将参考基因组进行索引构建，生成多个辅助文件，其中常见的文件格式包括：

.fasta.bwt: BWT 索引文件，存储了参考基因组的转换和排序信息，用于快速比对和查找匹配位置。（基于BWT算法）
.fasta.sa: 后缀数组文件，用于快速搜索和定位子序列的位置。
.fasta.amb: 参考基因组的字母和碱基的频率信息，用于压缩和优化索引。

对sam文件进行操作，并call SNP生成.vcf文件

#筛选比对上的高质量reads
samtools view -F4 -q1 -b map.sam -o map.bam#将低质量值去掉，并将sam转换为二进制的bam以方便计算
samtools sort -@ 2 map.bam map.sortP#使用两个线程来对bam进行排序（字典序），以方便后续计算
samtools index map.sort.bam#构建索引方便计算
#统计比对reads数
samtools view -c map.sortP.bam
#统计未比对上的reads数
samtools view -c -f 4 map.sam
#统计比对到正链上的reads数
samtools view -c -F 16 map.sam
#获取properly-paired的reads数
samtools view -f2 -F 256 -c map.sortP.bam
#查看每个位置碱基比对或错配的情况
samtools mplieup -f ref.fa -Q 13 map.sortP.bam | less
#(Q表示过滤掉低质量的测序碱基)

#call SNP(在上面结果没有问题的情况下，)
samtools mpileup -f ref.fa -Q 13 -vu map.sortP.bam > snp.vcf

以上是对一组双末端测序数据的call SNP操作，如果要进行正式的分析，我们往往需要懂得如何批量处理数据。我们可以将一组数据的分析流程与工具明确，将参数设置成自动，并使用循环或并行计算进行全部操作：

#!/bin/bash
data_dir="存放数据的绝对路径"
for group in group1 group2 group3
do
    # fastp 质控和过滤
    fastp -i "$data_dir/${group}_R1.fastq" -I "$data_dir/${group}_R2.fastq" \
    -o "$data_dir/${group}_filtered_R1.fastq" -O "$data_dir/${group}_filtered_R2.fastq"

    # bwa 进行比对
    bwa mem reference.fasta "$data_dir/${group}_filtered_R1.fastq" "$data_dir/${group}_filtered_R2.fastq" | samtools sort -o "$data_dir/${group}.sorted.bam"

    # 利用 samtools 做 SNP calling
    samtools mpileup -uf reference.fasta "$data_dir/${group}.sorted.bam" | bcftools call -mv > "$data_dir/${group}.vcf"
done