SeqKit教程：生物序列处理与分析实战指南
项目地址: https://gitcode.com/gh_mirrors/se/seqkit 实时监控病毒测序比对结果
在时间紧迫的疫情监测场景中，SeqKit提供了一系列实时处理FASTQ/FASTA和BAM格式数据的子命令，包括watch、fish、scat和bam等。这些功能为构建实时分析管道提供了强大支持。
下面我们展示一个完整的分析管道示例，该管道在下载病毒扩增子测序数据的同时，使用minimap2将其比对到参考序列上：
#!/bin/bash
set -o nounset
# 定义参考序列和数据URL
REF_URL="https://www.ncbi.nlm.nih.gov/sviewer/viewer.cgi?tool=portal&save=file&log$=seqview&db=nuccore&report=fasta&id=1798174254&extrafeat=null&conwithfeat=on&hide-cdd=on&ncbi_phid=CE8B108356DDCF110000000005B10489"
DATA_URL="http://ftp.sra.ebi.ac.uk/vol1/fastq/SRR145/055/SRR14560555/SRR14560555_1.fastq.gz"
# 构建实时分析管道
minimap2 -K 10M -t 8 -ax map-ont \
<(wget -q -O - "$REF_URL") \
<(wget -q -O - "$DATA_URL") 2>/dev/null \
| samtools view -F 2304 -q 1 -b - \
| seqkit bam -O aln_len.pdf -f RefAln -m 310 -M 410 -p 5000 -x - \
| seqkit bam -O aln_acc.pdf -f Acc -p 10000 -x - \
| samtools sort -o viral_artic_sorted.bam -
# 显示最终比对统计信息
seqkit bam -s -k viral_artic_sorted.bam
该管道实现了以下功能：
- 使用
minimap2进行实时比对 - 只保留主比对且比对质量大于1的记录
- 筛选比对参考长度在310-420之间的比对结果
- 每5000条记录显示一次比对长度分布直方图并保存为PDF
- 每10000条记录显示一次比对准确率分布直方图
- 最终输出排序后的BAM文件
序列去重与嵌套序列处理
在处理基因组组装结果时，经常需要去除重复序列和嵌套序列。以下是一个完整的处理流程：
示例数据
>big
ACTGACGATCGATACGCAGCACAGCAG
>small_in_big_rc
TGCTGCGTATCG
>small2
ACTACGACTACGACT
>small2_alias
ACTACGACTACGACT
>small2_rc
AGTCGTAGTCGTAGT
>another
ACTAACGA
处理步骤
- 去除完全相同的序列：
seqkit rmdup -s -i contigs.fa -o contigs.uniq1.fa
- 识别嵌套序列：
seqkit locate -M -f contigs.uniq1.fa contigs.uniq1.fa -o match.tsv
- 提取嵌套序列ID：
sed 1d match.tsv | awk '$2 != $1' | cut -f 2 | tee nested.txt
- 去除嵌套序列：
seqkit grep -v -f nested.txt contigs.uniq1.fa -o contigs.uniq2.fa
最终结果将只保留非重复且非嵌套的序列。
大型基因组处理技巧
处理人类基因组等大型序列时，SeqKit提供了多种优化方法：
1. 构建FASTA索引
seqkit faidx --id-regexp "^(.+)$" hsa.fa -o hsa.fa.seqkit.fai
2. 按序列长度排序
seqkit sort --by-length --reverse --two-pass hsa.fa > hsa.sorted.fa
3. 随机打乱序列顺序
seqkit shuffle hsa.fa --two-pass > hsa.shuffled.fa
4. 按染色体提取子序列
seqkit subseq -r 1:10 --chr X --chr Y hsa.fa
5. 根据GTF文件提取CDS序列
seqkit subseq --gtf Homo_sapiens.GRCh38.84.gtf.gz --chr X --feature cds hsa.fa > chrX.gtf.cds.fa
污染序列去除实战
- 首先获取污染序列的ID列表
- 使用
grep命令去除污染序列：
seqkit grep -f contaminate.list -v reads_1.fq.gz -o reads_1.clean.fq.gz
seqkit grep -f contaminate.list -v reads_2.fq.gz -o reads_2.clean.fq.gz
miRNA发夹结构分析
数据准备
从miRBase数据库下载hairpin.fa.gz文件
基本统计
seqkit stat hairpin.fa.gz
重复发夹结构分析
seqkit rmdup -s -i hairpin.fa.gz -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt
按物种分类
seqkit split hairpin.fa.gz -i --id-regexp "^([\w]+)\-" --two-pass
人类miRNA长度分布分析
seqkit grep -r -p '^hsa' hairpin.fa.gz | seqkit fx2tab -l | cut -f 4 | csvtk -H plot hist --xlab Length --title "Human pre-miRNA length distribution"
通过本教程，您已经掌握了SeqKit在生物序列处理中的多种实用技巧，从实时数据分析到大型基因组处理，再到miRNA特征分析。这些方法可以灵活应用于各种生物信息学分析场景。
创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
