# 實做測試
# **一.在linux系統下設置個人環境。**

1.我曾在几個月前因為對Linux充滿興趣，所以就提前在自己的win10系統上下載VMware虛擬機，然後并且從中下載Ubuntu64操作系統。
具體的方法我參考了 https://b23.tv/3838H1K 中的教學視頻。
主要是下載VMware軟體以及下載Ubuntu64系統的鏡像文件。

2.除此之外，我還通過查閲資料，創建了能成功連接win10系統和虛擬機Ubuntu系統的共享文件夾。
具體參考了 https://b23.tv/w8nzs2i 中的教學視頻。這個共享文件夾非常有用，可以讓我在win10和ubuntu中實時傳遞文件，這是雙系統所不具備的優點。

note：目前在使用VMware軟體中我曾經遇到了重新下載失敗的情況，後來查閲資料瞭解到，VMware卸載不幹淨，從而導致安裝失敗。而我在網路上找了很多方法，都覺得太過繁瑣。最後決定直接將系統重裝（格式化系統盤的所有數據，所幸我的數據都沒有放在系統磁盤）。

# **二.完成實作課程 https://datacarpentry.org/wrangling-genomics/ 中的內容**
#
EPISODES
Summary and Schedule
1. Background and Metadata
2. Assessing Read Quality
3. Trimming and Filtering
4. Variant Calling Workflow
5. Automating a Variant Calling

**1.Background and Metadata**
這個課程背景介紹了一個長期演化實驗，使用大腸桿菌（E. coli）菌株進行研究。該實驗在葡萄糖有限的基礎培養液中培養大腸桿菌族群超過50,000個世代。在實驗的不同時間點（5000、15000和50000個世代），取得了三個樣本。該實驗的目的是探索菌群的適應性變化，特別是涉及到能夠利用檸檬酸鹽的自發突變體（Cit+）和超突變性的變化

**2.Assessing Read Quality**
這個步驟進展非常順利，目前已經順利完成。過程中曾出現FastQC并未成功安装的問題。这个问题主要是因为在安装过程中有几个软件包无法下载。後來查閲資料尝试运行 sudo apt-get update 命令来更新我的软件源列表，然后再次尝试安装 FastQC后隨即成功。

note：我也通過查閲資料瞭解到質量得分（Quality Score）是在基因定序（DNA Sequencing）中非常重要的一個概念。在FASTQ文件中，每一個讀取（Read，即一段序列數據）都會有對應的質量得分。每一個質量得分對應一個讀取中的一個核苷酸，代表這個核苷酸的定序誤差概率。質量得分高，代表這個核苷酸的定序誤差概率低，反之亦然。

同時也瞭解到FastQC是一款生物信息學工具，主要用於品質控制(QC)的高通量序列數據分析。它可以從每個序列文件讀取數據，並產生一個由多個模塊組成的品質控制報告，每個模塊都可以幫助識別數據中可能存在的不同類型的問題。

### **3.Trimming and Filtering**
這個步驟目前順利完成。過程中也出現了跟上個步驟同樣的問題trimmomatic為成功安裝，但後續都有得到解決。以下是運行成功的截圖

大致的流程總結是：
1.下載並安裝了Trimmomatic，**並將其添加到環境變量中。**（添加到環境變量中非常重要，會影響後續的步驟）
2.在進行修剪和過濾之前，提前準備好FASTQ格式的原始讀取文件。
3.使用Trimmomatic命令來修剪和過濾我的讀取。需要指定一系列的參數和選項，包括讀取的輸入和輸出文件，修剪和過濾的標準，以及適配器序列等
4.Trimmomatic在結果中提供了詳細的修剪和過濾過程信息，包括讀取的總數，以及保留下來的讀取數。

note：
大概瞭解了有關質量控制（quality control）這個概念。
同時也大概理解了有關雙向讀取，前向讀取，反向讀取等概念。

**4.Variant Calling Workflow**


**problem：**

在上面的步驟中我遇到了困難。首先是解碼失敗。後來通過sub.tar.gz指令查詢到sub.tar.gz实际上是一个ASCII文本的HTML文档，而不是一个有效的tar或gzip压缩文件。这就解释了为什么我在尝试使用tar命令解压缩该文件时遇到了错误。
後來我嘗試打開Figshare 的网站并且搜索14418248文件的唯一的ID，但還是失敗了。猜測這個鏈接文件未公開或者已經被刪除。

**solution**
後來我認真閲讀題目瞭解到这个步骤在实际的分析中是非必需的，由于这些文件小于原始的FASTQ文件，所以处理它们的速度更快。这使得我能更快地运行和测试的变异检测工作流程，并看到结果。然而，在实际的分析中，我会使用原始的、完整的FASTQ文件，而不是这些小的子集。
因爲存在這個問題，所以我采用的方法是將上一個步驟中Trimming and Filtering產出的fastq文件選擇兩個SRR2584863_1.trim.fastq.gz 和 SRR2584863_2.trim.fastq.gz文件並將其複製到新創建的文件夾trimmed_fastq_small目錄下從而替代從網站下載的FASTQ的子集文件，從而能夠繼續以下的步驟。

**Variant Calling Workflow這一步驟相對於前面2，3兩個步驟來説要複雜得多。
所以要大概梳理一下其中的流程：**

**1.Alignment to a reference genome**
a.Indexing the reference genome
b.Aligning the reads to the reference genome

**2.alignemet cleanup**
a.SAM convert to BAM format
b.Sort BAM file by coordinates

**3.Variant calling**
a.Calculate the read coverage of positions in the genome
b.Detect the single nucleotide variants (SNVs)
c.Filter and report the SNV variants in variant calling format (VCF)

在經過上面一系列複雜的流程后，最後的呈現結果會是這樣：

note：用對齊工具（如BWA-MEM算法）將序列數據（讀取）對齊到參考基因組。這個步驟會生成一個包含對齊信息的SAM或BAM文件。
$ bwa mem ref_genome.fasta input_file_R1.fastq input_file_R2.fastq >output.sam 在這個指令中，需要注意的是要注意讀取的方向和配對信息。如果你的數據是配對端的（也就是說，每個讀取都有一個配對的讀取）。

note：在這裏我大致瞭解到了SAM與DAM這兩種資料形態的特點。BAM 文件是 SAM 的二進制版本并且讀寫速度更改，可以被索引。 所以我們需要將SAM 文件轉換為 BAM 文件進行存儲和後續分析 轉換的工具是通過****samtools****來實現的。

# **5.Automating a Variant Calling**（重點步驟）
這一步驟主要就是將前面所有的步驟統整結合起來，并且編寫兩個脚本來自動化地完成上述的步驟。

首先兩個脚本分別是
1.Analyzing quality with FastQC
2.Automating the rest of our variant calling workflow

我將上述的脚本厘清并且做好了注釋
----------------------------------------------
如果有错误发生，退出脚本
set -e

切换到 untrimmed_fastq/ 目录
cd ~/dc_workshop/data/untrimmed_fastq/

输出状态消息，并在所有.fastq文件上运行FastQC
echo "Running FastQC ..."
fastqc *.fastq*

创建目录来存放FastQC的结果
mkdir -p ~/dc_workshop/results/fastqc_untrimmed_reads

输出状态消息，并移动所有.zip和.html文件到结果目录
echo "Saving FastQC results..."
mv *.zip ~/dc_workshop/results/fastqc_untrimmed_reads/
mv *.html ~/dc_workshop/results/fastqc_untrimmed_reads/

切换到存储结果的目录
cd ~/dc_workshop/results/fastqc_untrimmed_reads/

输出状态消息，并遍历所有.zip文件进行解压
echo "Unzipping..."
for filename in *.zip
do
unzip $filename
done

输出状态消息，并将所有摘要文件连接成一个输出文件
echo "Saving summary..."
cat */summary.txt > ~/dc_workshop/docs/fastqc_summaries.txt

---------------------------------------------------------------

花時間去整理瞭解后，其實感覺這脚本並不複雜，就是一段簡單的程式碼。
不過這裏需要注意的是这个脚本中並没有明确的输入文件，但它依赖于操作系统的当前目录以及这个目录下特定扩展名（如.fastq、.zip和.html）的文件。具体来说，这个脚本做了以下事情：
*使用通配符.fastq来指定所有以.fastq为扩展名的文件作为FastQC的输入。
使用通配符.zip和.html来移动所有.zip和.html文件到指定的目录。
使用通配符.zip来找出所有的.zip文件并解压。
使用通配符/summary.txt来找出所有的summary.txt文件并连接

**至於2.Automating the rest of our variant calling workflow這個脚本要更加複雜得多**
variant calling workflow has the following steps:

1.Index the reference genome for use by bwa and samtools.
2.Align reads to reference genome.
3.Convert the format of the alignment to sorted BAM, with some intermediate steps.
4.Calculate the read coverage of positions in the genome.
5.Detect the single nucleotide variants (SNVs).
6.Filter and report the SNVs in VCF (variant calling format).

----------------------------------------------------------------------------

# 设置脚本在遇到错误时退出
set -e

# 切换工作目录到~/dc_workshop/results
cd ~/dc_workshop/results

# 将参考基因组的文件路径存储在变量genome中
genome=~/dc_workshop/data/ref_genome/ecoli_rel606.fasta

# 使用BWA软件对参考基因组进行索引
bwa index $genome

# 创建用于存储不同类型结果文件的目录
mkdir -p sam bam bcf vcf

# 对~/dc_workshop/data/trimmed_fastq_small/*_1.trim.sub.fastq路径下的所有fastq文件进行处理
for fq1 in ~/dc_workshop/data/trimmed_fastq_small/*_1.trim.sub.fastq
do
# 显示正在处理的文件名
echo "working with file $fq1"

# 提取fastq文件的基本名，存储在变量base中
base=$(basename $fq1 _1.trim.sub.fastq)
echo "base name is $base"

# 定义输入和输出文件的路径
fq1=~/dc_workshop/data/trimmed_fastq_small/${base}_1.trim.sub.fastq
fq2=~/dc_workshop/data/trimmed_fastq_small/${base}_2.trim.sub.fastq
sam=~/dc_workshop/results/sam/${base}.aligned.sam
bam=~/dc_workshop/results/bam/${base}.aligned.bam
sorted_bam=~/dc_workshop/results/bam/${base}.aligned.sorted.bam
raw_bcf=~/dc_workshop/results/bcf/${base}_raw.bcf
variants=~/dc_workshop/results/vcf/${base}_variants.vcf
final_variants=~/dc_workshop/results/vcf/${base}_final_variants.vcf

# 使用BWA软件进行序列比对，生成.sam格式的比对结果
bwa mem $genome $fq1 $fq2 > $sam

# 使用samtools将.sam格式的文件转换为.bam格式
samtools view -S -b $sam > $bam

# 对.bam格式的文件进行排序
samtools sort -o $sorted_bam $bam

# 对排序后的.bam文件进行索引
samtools index $sorted_bam

# 使用bcftools计算每个位置的读取覆盖度
bcftools mpileup -O b -o $raw_bcf -f $genome $sorted_bam

# 使用bcftools进行变异检测，生成.vcf文件
bcftools call --ploidy 1 -m -v -o $variants $raw_bcf

# 使用vcfutils进行变异过滤，生成最终的.vcf文件
vcfutils.pl varFilter $variants > $final_variants
done

最後是結果的呈現：

**question：**
在最後結果中，我後來再次檢查確認，發現run_variant_calling.sh這個脚本并沒有按照我所期望的那樣能夠自動化地產出vcf檔案。後來我仔細研究網站上的脚本意識到一個問題：

**season：**
在第二个脚本 run_variant_calling.sh 中，脚本期望的输入文件是已经被处理过的 fastq 文件（命名模式为 ${base}_1.trim.sub.fastq 和${base}_2.trim.sub.fastq），而不是原始的 fastq 文件。这个“处理过”的文件应该是对原始 fastq 文件进行質量控制后的结果（例如，修剪、滤波等）。
如果我的第一个脚本没有生成这样的处理过的 fastq 文件，那么你的第二个脚本run_variant_calling.sh 可能就无法找到它需要的输入文件。这就是我在运行第二个脚本时遇到问题的原因。

**solution：**

面對這個問題：有兩個主要的解決方法，我需要找到或生成适当的输入文件供第二个脚本使用。这可能意味着我需要修改網站上的第一个脚本，让它不仅进行质控，还要进行一些处理步骤，并生成新的 fastq 文件。或者，你可能需要找到一个额外的步骤或脚本来进行这个处理。
另一种解决办法是，如果的数据已经是足够高质量的，或者不需要进行任何预处理步骤，我可以修改的第二个脚本，让它接受原始的 fastq 文件作为输入。这就意味着我需要修改run_variant_calling.sh这部分的代码：
目前我的想法是直接采用第一種方法。

**於是我又重新寫了一個名叫process_qc.sh的脚本使其能夠通過 Trimmomatic程序對其進行質量控制，從而刪除質量不好的序列，并且輸出的文件能夠作爲本run_variant_calling.sh的輸入文件。**

#!/bin/bash

# 设置脚本在遇到错误时退出
set -e

# 切换工作目录到~/dc_workshop/data/untrimmed_fastq
cd ~/dc_workshop/data/untrimmed_fastq

# 对 ~/dc_workshop/data/untrimmed_fastq/ 路径下的所有 fastq 文件进行处理
for fq in *_1.fastq.gz
do
# 打印正在处理的文件名
echo "Processing file $fq"

# 提取 fastq 文件的基本名，存储在变量 base 中
base=$(basename $fq _1.fastq.gz)
echo "base name is $base"

# 定义输入和输出文件的路径
fq1=${base}_1.fastq.gz
fq2=${base}_2.fastq.gz
outpath=~/dc_workshop/data/trimmed_fastq_small
trimmed_fq1=${outpath}/${base}_1.trim.sub.fastq
trimmed_fq2=${outpath}/${base}_2.trim.sub.fastq
unpaired_fq1=${outpath}/${base}_1un.trim.sub.fastq
unpaired_fq2=${outpath}/${base}_2un.trim.sub.fastq

# 使用 Trimmomatic 对 fastq 文件进行质控，并生成处理后的 fastq 文件
TrimmomaticPE \
-threads 1 \
-phred33 \
$fq1 $fq2 \
$trimmed_fq1 $unpaired_fq1 \
$trimmed_fq2 $unpaired_fq2 \
ILLUMINACLIP:/usr/share/trimmomatic/NexteraPE-PE.fa:2:30:10 \
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
done

并且這個脚本修改了原先read_qc.sh這個脚本所沒能夠解決的問題。 process_qc.sh對原始的fastq檔案進行了處理并且生成的輸出文件可作爲脚本run_variant_calling.sh的輸入文件。進而能夠全自動的完成上述步驟。

經過大量的時間和精力調試，process_qc.這個脚本終於能夠成功運行，下圖為proces_qc脚本的執行結果。

可以看到process_qc脚本成功運行并且成功產出了fastq文件，能夠為後續的脚步運行提供輸入文件。

緊接著我就開始執行run_variant_calling.sh這個脚本文件。大概運行了20分鐘左右，最後成功執行完畢。并且成功產生了三個vcf文件，分別為SRR2589044，SRR2584863， SRR2584866，分別對應這E.coil對應的第5000，15000，50000 generation。

-------------------------------------------------------------------------


於是緊接著我通過簡單的工具 **bcftools **去計算這三個vcf文件中每個染色體上的變異數量。最後得出了以下的結果。

可以看出E.coil與5000，15000，50000 generation對應的變異數量分別是11，29，794個變異數量。可以看出隨著時間的推移，E.coil的變異數量不斷增大。

--------------------------------------------------------------
至於爲什麽與上面網路上的教學結果不太一樣，我猜測是由於我是用trimmed工具進行quality control操作時，調用的參數不太一樣，從而導致最後產生的變異數量不一樣。

看到上面的結果讓我非常開心，因爲這算是成功完成上面的實作測試，成功通過自動化脚本產出三個vcf文件，并且分析出E.coil染色體上的變異數量。但是這並不能全面的反映出E.coil的變異程度。通過查閲資料我瞭解到Python 的 pysam 库就可以用来读取和处理 VCF 文件，從而讀取更加詳細複雜的訊息。

# 三.通過python解析VCF檔案

我用python并且下載pysam套件，成功讀取出vcf檔案，并且將其轉換成pandas的dataframe格式。

單單看變異數量我還是不夠全面，我想進一步查看這些變異對應的染色體位置分佈。於是我想到了或許可以通過用python將E.coil的基因組位置坐標平均分成100份，每一份統計對應的變異次數，從而構建一個10×10的二維熱圖。我從大致瞭解到E.coil genome總共有463萬個基因堿基，所以每一份對應4.63萬個堿基坐標。

下面是我的程式碼

并且我成功繪製出了對應5000，15000，50000generation對應的熱圖




除了查看變異分佈之外，我還想查看變異類型。於是我根據REF（參考序列）和 ALT（變異序列）這兩個column可以將變異類型分成3種，分別是SNPs，MNVs，Indels.并且繪製餅圖。以下是繪製的結果。


-----------------------------------------------

-------------------------------------------------

-------------------------------------------------------

除了查看變異分佈之外，我還通過繪製boxplot去繪製出變異質量的分佈。
以下是程式碼以及對應的結果：



# 改進：通過python脚本一鍵視覺化解析vcf檔案

做好了上面一系列的分析之後，我突然想到，是否也可以像跟上面的脚本一樣。只要執行python檔案就可以一鍵輸出上面所有的數據分析圖片，并且將所有的圖片放入同一個文件夾内。於是照著這個思路我將python code進一步進行了改進。以下是我成功編寫的python脚本

import glob
import pysam
import numpy as np
import pandas as pd
import os

####################################
#vcf convert to dataframe
###################################

def create_dataframe_from_vcf(file_path):
# 讀取 VCF 文件
vcf_reader = pysam.VariantFile(file_path)

# 建立一個字典來存儲信息
records_dict = {'CHROM': [], 'POS': [], 'REF': [], 'ALT': [], 'QUAL': [] ,'INFO': []}

# 遍歷 VCF 文件的每一個變異
for record in vcf_reader:
records_dict['CHROM'].append(record.chrom)
records_dict['POS'].append(record.pos)
records_dict['REF'].append(record.ref)
records_dict['ALT'].append(record.alts)
records_dict['QUAL'].append(record.qual)
records_dict['INFO'].append(record.info)

# 將字典轉換為 pandas DataFrame
df = pd.DataFrame(records_dict)

return df

import glob

# 定义文件路径模式
vcf_file_pattern = "/home/clarence/dc_workshop/results/vcf/*_final_variants.vcf"

# 使用 glob 查找所有匹配的文件
vcf_files = glob.glob(vcf_file_pattern)

# 创建一个字典，键是文件名，值是对应的 DataFrame
df_dict = {}
for vcf_file in vcf_files:
# 提取文件名作为 DataFrame 的键
file_key = vcf_file.split('/')[-1].split('_final')[0]

# 创建 DataFrame 并存储在字典中
df_dict[file_key] = create_dataframe_from_vcf(vcf_file)

# 创建data_analysis_report文件夹
folder_name = "data_analysis_visualization"
if not os.path.exists(folder_name):
os.makedirs(folder_name)

######################################
#Variant_Distribution_Heatmap
#####################################

import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import os

def create_heatmap(df, title):
num_bins=100
total_length=4629225
bin_size=total_length//num_bins
mutation_counts = [0]*num_bins

for index, row in df.iterrows():
bin_index = min(row['POS'] // bin_size, num_bins-1)
mutation_counts[bin_index] += 1

mutation_counts_matrix = np.array(mutation_counts).reshape((10,10))

sns.heatmap(mutation_counts_matrix, annot=True, cmap='YlGnBu')
plt.title(title)

# 在保存图像时指定文件夹
plt.savefig(os.path.join(folder_name, title + "_Variant_Distribution_Heatmap.png"))
plt.close()

for title, df in df_dict.items():
create_heatmap(df, title)

import matplotlib.pyplot as plt
import os

# 定义一个函数来判断变异类型
def classify_variant(row):
if 'SVTYPE' in row['INFO']:
return 'SV'
if len(row['REF']) == 1 and len(row['ALT'][0]) == 1:
return 'SNP'
elif len(row['REF']) > 1 and len(row['ALT'][0]) > 1:
return 'MNV'
else:
return 'Indel'

#################################
#Variant_Type_PieChart
################################

def create_pie(df, title):
# 创建新的列 "Mutation Type"
df['Mutation Type'] = df.apply(classify_variant, axis=1)

# 计算每种变异类型的数量
mutation_counts = df['Mutation Type'].value_counts()

# 画饼图
plt.figure(figsize=(8, 8))
plt.pie(mutation_counts, labels=mutation_counts.index, autopct='%1.1f%%')
plt.title(title + ' Variation Types')

# 在保存图像时指定文件夹
plt.savefig(os.path.join(folder_name, title + "_Variant_Type_PieChart.png"))
plt.close()

# 假设 df_dict 是一个包含多个 dataframe 的字典，
# 其中每个 dataframe 都有一个相应的标题（用作文件名）
for title, df in df_dict.items():
create_pie(df, title)

###########################
#Variant_Quality_Boxplot
###########################

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

# 假设df_dict是一个包含多个dataframe的字典，键是每个dataframe的来源
for source, df in df_dict.items():
df['source'] = source

# 将所有dataframe合并为一个
combined_df = pd.concat(df_dict.values())

# 创建boxplot
sns.set_style("whitegrid")

plt.figure(figsize=(10,8))
sns.boxplot(x='source', y='QUAL', data=combined_df, width=0.4)

plt.title('Variant Quality Boxplot')
plt.xlabel('Group')
plt.ylabel('Quality')

# Save the figure before showing it
plt.savefig(os.path.join(folder_name, 'Variant_Quality_Boxplot.png'))

# Close the figure
plt.close()

通過以上這個名叫data_visual_analysis.py的脚本，使我能夠一鍵解析該路徑下所有的vcf檔案並將其繪製成圖片從而視覺化分析,并且在相同路徑下創建一個名叫**data_analysis_visualization**的文件夾，將所有的分析圖片放入該文件夾内。如下圖所示：

并且測試完脚本沒問題后，我可以直接在linux的終端去運行python脚本，非常地方便。

紅色方框内的圖片是 vcf_visual_analysis.py脚本產生的針對相同路徑下的vcf檔案解析并且產生的圖片。圖片的内容就是上面在vscode數據分析時呈現的heatmap，piechart，boxplot等結果。

----------------------------------------------------------

# 脚本流程優化

所以流程全部流程下來總共有三個脚本，分別是process_qc.sh,variant_calling.sh以及vcf_visual_analysis.py的python脚本。

## **我將上述所有的三個脚本放入了Google drive雲端,以下為鏈接：**

## https://drive.google.com/drive/folders/12y-rhVOUJkq7MODZsw5-SjOJsjAxMTEd?usp=drive_link
通過對於流程的分析，後來我認爲可以再對process_qc.sh這個脚本進行優化。
原先的process_qc.sh脚本主要是進行 fastq read和trimmed步驟。
後來我仔細思考了一下，或許可以對比fasq文件通過trimmed進行quality control前後都進行fastq_read步驟進而可以分析前後對比的結果。

於是process_qc的流程由原來的
fastq_read → quality control
改進成
fastq_read → quality control → fastq_read
所以修改后的process_qc脚本代碼如下：

#!/bin/bash

# 创建需要的目录
mkdir -p ~/dc_workshop/results/fastqc_untrimmed_reads
mkdir -p ~/dc_workshop/results/fastqc_trimmed_reads

# 设置脚本在遇到错误时退出
set -e

# 切换工作目录到~/dc_workshop/data/untrimmed_fastq
cd ~/dc_workshop/data/untrimmed_fastq

# 对 ~/dc_workshop/data/untrimmed_fastq/ 路径下的所有 fastq 文件进行处理
for fq in *_1.fastq.gz
do
# 打印正在处理的文件名
echo "Processing file $fq"

# 提取 fastq 文件的基本名，存储在变量 base 中
base=$(basename $fq _1.fastq.gz)
echo "base name is $base"

# 定义输入和输出文件的路径
fq1=${base}_1.fastq.gz
fq2=${base}_2.fastq.gz
outpath=~/dc_workshop/data/trimmed_fastq_small
trimmed_fq1=${outpath}/${base}_1.trim.sub.fastq
trimmed_fq2=${outpath}/${base}_2.trim.sub.fastq
unpaired_fq1=${outpath}/${base}_1un.trim.sub.fastq
unpaired_fq2=${outpath}/${base}_2un.trim.sub.fastq

# 运行 FastQC 对原始数据进行质量分析，生成 HTML 报告
fastqc -o ~/dc_workshop/results/fastqc_untrimmed_reads $fq1 $fq2

# 使用 Trimmomatic 对 fastq 文件进行质控，并生成处理后的 fastq 文件
TrimmomaticPE \
-threads 1 \
-phred33 \
$fq1 $fq2 \
$trimmed_fq1 $unpaired_fq1 \
$trimmed_fq2 $unpaired_fq2 \
ILLUMINACLIP:/usr/share/trimmomatic/NexteraPE-PE.fa:2:30:10 \
LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

# 运行 FastQC 对处理后的数据进行质量分析，生成 HTML 报告
fastqc -o ~/dc_workshop/results/fastqc_trimmed_reads $trimmed_fq1 $trimmed_fq2
done

爲了更好地理解脚本的流程，我畫了一張對應的流程圖，裏面有脚本的輸入文件和輸出文件，以及使用的工具。這張圖可以比較直觀地看出脚本的流程。
簡單將就是三個步驟：fastq_read → quality control → fastq_read

------------------------------------------------

緊接著我又嘗試做出了variant_calling.sh的流程圖。相對而言variant_calling.脚本的流程就要複雜許多。

以下是詳細步驟：

1.BWA索引: 使用BWA对参考基因组进行索引。

2.BWA序列比对: 使用BWA将FASTQ文件与参考基因组进行比对，产生SAM文件。

3.SAM到BAM的转换: 使用samtools将SAM文件转换为BAM文件。

4.BAM文件排序和索引: 对BAM文件进行排序并创建索引。

5.读取覆盖度计算: 使用bcftools计算每个位置的读取覆盖度，输出BCF文件。

6.变异检测: 使用bcftools进行变异检测，将BCF文件转换为VCF文件。

7.变异过滤: 使用vcfutils进行变异过滤，生成最终的VCF文件。