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)

————————————————————————————————————

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 219,366评论 6 508
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 93,521评论 3 395
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 165,689评论 0 356
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,925评论 1 295
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,942评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,727评论 1 305
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,447评论 3 420
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,349评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,820评论 1 317
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,990评论 3 337
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,127评论 1 351
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,812评论 5 346
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,471评论 3 331
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,017评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,142评论 1 272
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,388评论 3 373
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 45,066评论 2 355

推荐阅读更多精彩内容