shell脚本统计reads数量和碱基数量
1.github下载命令
git clone 【网页】
2.安装miniconda
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-py38_4.9.2-Linux-x86_64.sh
bash Miniconda3-py38_4.9.2-Linux-x86_64.sh
3.查看gzip压缩文件中的前几行
zcat filename.gz | head -n 10
4.查看gzip压缩文件中的某个字符的个数
zcat yourfile.gz | grep -o '@' | wc -l
5.seqkit
seqkit是一个用于处理生物信息学序列数据的跨平台工具包，它支持多种与 FASTA/Q 文件相关的操作。以下是一些常用的seqkit命令及其用法：
查看帮助信息：
seqkit
查看版本：
seqkit --version
统计序列数量：
seqkit stat in.fasta
查看一个fastq文件中的每条序列的长度和CG值
seqkit fx2tab -l -g -n -i -H in.fasta
提取序列子集：
seqkit sample in.fasta -n 10 > out.fasta
序列过滤：
根据序列长度过滤：
seqkit filter in.fasta --min-len 100 -o out.fasta
根据序列ID过滤：
seqkit grep -p "序列ID模式" in.fasta > out.fasta
序列转换：
DNA 转 RNA：
seqkit seq -d2r in.fasta > out.fasta
RNA 转 DNA：
seqkit seq -r2d in.fasta > out.fasta
序列补码：
seqkit seq -p in.fasta > out.fasta
序列反转：
seqkit seq -r in.fasta > out.fasta
序列去gap：
seqkit seq -g in.fasta > out.fasta
序列格式转换：
将 FASTQ 转换为 FASTA：
seqkit seq in.fastq > out.fasta
将压缩的 FASTA/Q 文件转换为非压缩格式：
seqkit seq in.fasta.gz > out.fasta
序列提取：
提取具有特定特征的序列：
seqkit amplicon -f 1 -r 2 -s 100 in.fasta > out.fasta
提取 FASTQ 文件中的特定序列：
seqkit seq in.fastq -p -t DNA > out.fasta
序列比对：
seqkit map in.fasta -r ref.fasta > map.tsv
序列排序：
seqkit sort in.fasta > sorted.fasta
序列合并：
seqkit merge in.fasta1 in.fasta2 > merged.fasta
序列打乱：
seqkit shuffle in.fasta > shuffled.fasta
6.用seqkit来计数
seqkit工具可以用来查看 FASTA 或 FASTQ 文件中的序列数（即 reads）。在生物信息学中，一个 read 通常指的是测序得到的一个独立的 DNA 或 RNA 序列片段。在 FASTQ 文件中，每个 read 由四行组成：序列标识符、序列数据、+号（可选）和质量分数。
要使用seqkit查看 FASTQ 文件中的 reads 数量，可以使用stat命令，如下所示：
seqkit stat yourfile.fastq
这个命令将输出文件的一些统计信息，包括总的序列数（reads 数）。
如果你只想快速查看 reads 的数量，可以使用以下命令：
seqkit stat -f fastq yourfile.fastq | grep 'Total sequences'
这里-f fastq指定了输入文件是 FASTQ 格式的，然后通过管道将输出传递给grep命令，以筛选出包含“Total sequences”的行。
请注意，seqkit默认按照每四个行作为一个 read 来计数，因为 FASTQ 格式的每个 read 由四行组成。如果你的 FASTQ 文件被截断或者格式不正确，seqkit可能会报告一个异常的 reads 数量。如果需要处理压缩的 FASTQ 文件（如.fastq.gz或.fq.gz），seqkit会自动处理压缩，无需额外指定。
脚本
————————————————————————————————————
#!/bin/bash
# 定义输入文件夹和输出文件名
input_dir="/ifs1/01.RawData/03.Hic"
output_dir="/home/chensiyuan/out"
output_file="$output_dir/seqkit_stats_output.txt"
# 创建输出目录，如果不存在的话
mkdir -p "$output_dir"
# 清空输出文件，以防之前的内容影响
> "$output_file"
# 遍历指定文件夹下的所有.gz文件
for file in "$input_dir"/*.gz; do
# 提取文件名（不包含路径和扩展名）
filename=$(basename -- "$file")
filename_noext="${filename%.*}"
# 打印当前处理的文件名
echo "Processing file: $filename"
# 运行seqkit stats命令并将结果追加到输出文件中
seqkit stats "$file" >> "$output_file"
done
echo "All files processed. Output saved to $output_file"
————————————————————————————————————
import subprocess
import pandas as pd
import os
def get_seqkit_stats(file_path):
try:
result = subprocess.run(['seqkit', 'stat', file_path], capture_output=True, text=True, check=True)
return result.stdout
except subprocess.CalledProcessError as e:
print(f"An error occurred while processing {file_path}: {e}")
return None
def collect_stats(directory):
stats_list = []
for file_name in os.listdir(directory):
if file_name.endswith('.gz') :
file_path = os.path.join(directory, file_name)
stats = get_seqkit_stats(file_path)
if stats:
stats_list.append({
'File Name': file_name,
'Stats': stats
})
return stats_list
def save_to_excel(stats_list, output_file):
df = pd.DataFrame(stats_list)
df.to_excel(output_file, index=False)
directory_path = '/ifs1/01.RawData/01.HiFi'
output_excel_path = '/home/chensiyuan'
stats_list = collect_stats(directory_path)
save_to_excel(stats_list, output_excel_path)
print(f'All fasta/q files stats have been saved to {output_excel_path}')
————————————————————————————————————
import subprocess
import pandas as pd
import os
def get_seqkit_stats(file_path):
try:
result = subprocess.run(['seqkit', 'stat', file_path], capture_output=True, text=True, check=True)
return result.stdout.splitlines() # 修改：返回统计信息的列表
except subprocess.CalledProcessError as e:
print(f"An error occurred while processing {file_path}: {e}")
return None
def collect_and_save_stats(directory, output_file):
writer = pd.ExcelWriter(output_file, engine='openpyxl') # 创建ExcelWriter对象
stats_count = 0
for file_name in os.listdir(directory):
if file_name.endswith('.gz'):
file_path = os.path.join(directory, file_name)
stats = get_seqkit_stats(file_path)
if stats:
# 将统计信息转换为DataFrame
df = pd.DataFrame(stats_list, columns=['File Name', 'Stats'])
df.insert(0, 'File Name', file_name) # 在第一列插入文件名
df.to_excel(writer, sheet_name=file_name, index=False) # 写入Excel文件
stats_count += 1
print(f"Written stats for {file_name} to {output_file}")
writer.save() # 保存Excel文件
writer.close()
print(f'All stats have been saved to {output_file} with {stats_count} files.')
# 脚本主逻辑
directory_path = '/ifs1/01.RawData/01.HiFi'
output_excel_path = '/home/chensiyuan/stats_summary.xlsx' # 修改：指定完整的Excel文件名
collect_and_save_stats(directory_path, output_excel_path)
————————————————————————————————————
