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

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

SNP挖掘

质控

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进行双端测序样本修剪的示例命令，包含了一些常用的参数设置：

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

• <num_threads>：指定您要使用的线程数，根据您的机器配置和需要进行设置。

• <input1.fastq>和<input2.fastq>：替换为您的双端测序样本的输入文件路径。

• <output1_paired.fastq>和<output1_unpaired.fastq>：设置第一个样本的修剪结果输出路径。

• <output2_paired.fastq>和<output2_unpaired.fastq>：设置第二个样本的修剪结果输出路径。

• <adapter.fa>：如果适用，请提供您用于去除适配体序列的FASTA文件路径。如果不需要去除适配体序列，请删除整个ILLUMINACLIP:<adapter.fa>:2:30:10参数。

• 参数设置如LEADING:15和TRAILING:15表示移除序列开头和结尾低质量碱基的阈值，您可以根据实际情况调整这些值。

• SLIDINGWINDOW:4:15指定了滑动窗口的大小和最小平均质量要求，您也可以根据需要进行调整。

• MINLEN:36设置修剪后序列的最小长度，根据您的实验需求进行调整。

这只是一个基本的示例命令，您可以根据具体需求和数据情况进行参数调整和优化。确保在运行命令之前已经安装了Trimmomatic并正确设置了其环境变量。

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

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 是一个常用的基因组比对工具，它使用 Burrows-Wheeler Transform（BWT）索引算法来实现快速且高效的比对。

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

1. .fasta.bwt: BWT 索引文件，存储了参考基因组的转换和排序信息，用于快速比对和查找匹配位置。

2. .fasta.sa: 后缀数组文件，用于快速搜索和定位子序列的位置。

3. .fasta.amb: 参考基因组的字母和碱基的频率信息，用于压缩和优化索引。

第一步构建索引之后会生成.amb .pac .ann .bwt .sa等文件，具体情况如下：
.fasta.bwt文件:是基于BWT算法的压缩技术，基于BWT算法，可以利用其FL序列回溯构建出原参考基因组

.fasta.sa 文件（后缀数组）：后缀数组是一个有序的后缀索引列表，可以用于快速定位和搜索字符串中的子序列。利用 .fasta.sa 文件，可以实现以下操作：

• 子序列搜索：通过后缀数组可以快速找到特定子序列首字母的出现位置，然后通过 .fasta.bwt 文件进行进一步的匹配和定位。
• 出现次数统计：通过后缀数组，可以快速确定某个特定子序列在原始序列中出现的次数。

.fasta.pac 文件（前缀字典）：前缀字典是一个存储了原始序列数据的第一个字符的前缀信息的辅助文件。通过 .fasta.pac 文件，可以实现以下操作：

• 碱基匹配：利用 .fasta.pac 文件中的前缀信息，可以快速确定碱基对应的位置，帮助实现高效的碱基匹配操作，如 DNA 序列中的 A、C、G、T 碱基。
• 辅助索引：.fasta.pac 文件可以用作进一步优化索引的辅助数据结构，加速索引构建和索引查询的过程。

对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