seqkit+awk+sed+grep高级用法技巧合辑

生信与基因组学



软件安装与示例数据文件下载

ensembl最新版本下载网址： https://ftp.ensembl.org/pub/release-112/

# seqkit安装
conda install -c bioconda seqkit -y

# 格式转换
conda install csvtk -y

# 查看seqkit帮助
seqkit

1. seqkit 查看fa.gz和fq.gz序列文件

# 下载GRCh38 fasta文件
wget https://ftp.ensembl.org/pub/release-112/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
# 解压
# gzip -d Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# 下载GRCh38 gtf文件
wget https://ftp.ensembl.org/pub/release-112/gtf/homo_sapiens/Homo_sapiens.GRCh38.112.gtf.gz

seqkit可自动识别文件扩展名，而无需使用zcat查看.gz文件。

seqkit seq hairpin.fa.gz|less -S
# >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
# UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
# UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA

seqkit seq reads_1.fq.gz|less -S
# @HWI-D00523:240:HF3WGBCXX:1:1101:2574:2226 1:N:0:CTGTAG
# TGAGGAATATTGGTCAATGGGCGCGAGCCTGAACCAGCCAAGTAGCGTGAAGGATGACTGCCCTACGGG
# +
# HIHIIIIIHIIHGHHIHHIIIIIIIIIIIIIIIHHIIIIIHHIHIIIIIGIHIIIIHHHHHHGHIHIII


# 仅显示序列名称
seqkit seq Homo_sapiens.GRCh37.dna.primary_assembly.fa -n|head
# 1 dna:chromosome chromosome:GRCh37:1:1:249250621:1 REF
# 10 dna:chromosome chromosome:GRCh37:10:1:135534747:1 REF
# 11 dna:chromosome chromosome:GRCh37:11:1:135006516:1 REF

2. seqkit 压缩fasta文件

# -o 自动识别文件扩展名， 将fasta压缩为.gz格式
seqkit seq Homo_sapiens.GRCh37.dna.primary_assembly.fa.fasta \
-o Homo_sapiens.GRCh37.dna.primary_assembly.fa.fasta.gz

3. seqkit 统计fasta/fastq文件

# 单个文件统计
seqkit stats Homo_sapiens.GRCh37.dna.primary_assembly.fa
# file                                         format  type  num_seqs        sum_len  min_len       avg_len      max_len
# Homo_sapiens.GRCh37.dna.primary_assembly.fa  FASTA   DNA         84  3,101,804,739    4,262  36,926,246.9  249,250,621


# 统计全部fq.gz 和 fq.gz文件全部信息
seqkit stats *.f{a,q}.gz -a
# file               format  type  num_seqs    sum_len  min_len  avg_len  max_len   Q1   Q2   Q3  sum_gap  N50  N50_num  Q20(%)  Q30(%)  AvgQual  GC(%)
# reads_1.fq.gz      FASTQ   DNA      2,500    567,516      226      227      229  227  227  227        0  227        3   91.24   86.62    15.45  53.63
# reads_2.fq.gz      FASTQ   DNA      2,500    560,002      223      224      225  224  224  224        0  224        2   91.06   87.66    14.62  54.77


# 表格格式输出， 并通过csvtk转换为tab分割
seqkit stats *.f{a,q}.gz -T | csvtk pretty -t
# file                format   type   num_seqs   sum_len   min_len   avg_len   max_len
# -----------------   ------   ----   --------   -------   -------   -------   -------
# reads_1.fq.gz       FASTQ    DNA    2500       567516    226       227.0     229
# reads_2.fq.gz       FASTQ    DNA    2500       560002    223       224.0     225

4. seqkit 提取fasta指定染色体序列

# 提取1号染色体序列
seqkit grep -p 1 Homo_sapiens.GRCh37.dna.primary_assembly.fa \
-o Homo_sapiens.GRCh37.dna.primary_assembly.chr1.fa


# 查看>开头信息
cat Homo_sapiens.GRCh37.dna.primary_assembly.fa|grep '^>'
# 提取2号染色体10000-10050碱基序列
seqkit faidx Homo_sapiens.GRCh37.dna.primary_assembly.fa 2:10000-10050
#[INFO] create or read FASTA index ...
#[INFO] read FASTA index from Homo_sapiens.GRCh37.dna.primary_assembly.fa.fai
#[INFO]   84 records loaded from Homo_sapiens.GRCh37.dna.primary_assembly.fa.fai
#>2:10000-10050
#NCGTATCCCACACACCACACCCACACACCACACCCACACACACCCACACC

5. seqkit合并fastq文件夹全部fastq文件

seqkit scat -j 4 -f fastq_dir > all.fq

6. seqkit转换fasta为氨基酸序列文件

# 查看密码子表详细表, notincluding ambigugous codons
seqkit translate -l 1

# 查看密码子表详细表, including ambigugous codons
seqkit translate -L 1
# The Standard Code (transl_table=1)
# Source: https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes#SG1

# Initiation Codons:
#     ATG, CTG, TTG

# Stop Codons:
#     TAA, TAG, TGA

# Stranslate Table:
#     AAA: K, AAC: N, AAG: K, AAT: N
#     ACA: T, ACC: T, ACG: T, ACT: T
#     AGA: R, AGC: S, AGG: R, AGT: S
#     ATA: I, ATC: I, ATG: M, ATT: I

#     CAA: Q, CAC: H, CAG: Q, CAT: H
#     CCA: P, CCC: P, CCG: P, CCT: P
#     CGA: R, CGC: R, CGG: R, CGT: R
#     CTA: L, CTC: L, CTG: L, CTT: L

#     GAA: E, GAC: D, GAG: E, GAT: D
#     GCA: A, GCC: A, GCG: A, GCT: A
#     GGA: G, GGC: G, GGG: G, GGT: G
#     GTA: V, GTC: V, GTG: V, GTT: V

#     TAA: *, TAC: Y, TAG: *, TAT: Y
#     TCA: S, TCC: S, TCG: S, TCT: S
#     TGA: *, TGC: C, TGG: W, TGT: C
#     TTA: L, TTC: F, TTG: L, TTT: F

# --trim参数表示去除*
seqkit translate mouse-p53-cds.fna --trim
# >lcl|AB021961.1_cds_BAA82344.1_1 [gene=p53] [protein=P53] [protein_id=BAA82344.1] #[location=101..1273] [gbkey=CDS]
# MTAMEESQSDISLELPLSQETFSGLWKLLPPEDILPSPHCMDDLLLPQDVEEFFEGPSEA
# LRVSGAPAAQDPVTETPGPVAPAPATPWPLSSFVPSQKTYQGNYGFHLGFLQSGTAKSVM
# CTYSPPLNKLFCQLAKTCPVQLWVSATPPAGSRVRAMAIYKKSQHMTEVVRRCPHHERCS
# DGDGLAPPQHRIRVEGNLYPEYLEDRQTFRHSVVVPYEPPEAGSEYTTIHYKYMCNSSCM
# GGMNRRPILTIITLEDSSGNLLGRDSFEVRVCACPGRDRRTEEENFRKKEVLCPELPPGS
# AKRALPTCTSASPPQKKKPLDGEYFTLKIRGRKRFEMFRELNEALELKDAHATEESGDSR
# AHSSYLKTKKGQSTSRHKKTMVKKVGPDSD

7. 连接碱基序列

echo -e ">seq\nACGT-ACTGC-ACC"| seqkit seq -g -u
# >seq
# ACGTACTGCACC

8. 使用序列长度过滤fasta文件

-m 最小长度， -M 最大长度

cat hairpin.fa| seqkit seq -m 100 -M 1000| seqkit stats
# 文件格式类型num_seqs    sum_len min_len avg_len max_len
#-     FASTA   RNA     10,972 1560270      100    142.2      938

9. RNA和DNA序列的相互转换

# RNA to DNA, DNA to RNA 用--dna2rna参数
echo -e ">seq\nUCAUAUGCUUGUCUCAAAGAUUA" | seqkit seq --rna2dna
# >seq
# TCATATGCTTGTCTCAAAGATTA

10. 输出反向互补序列

seqkit seq hairpin.fa.gz -r -p
# >cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
# UCGAAGAGUUCUGUCUCCGGUAAGGUAGAAAAUUGCAUAGUUCACCGGUGGUAAUAUUCC
#AAACUAUACAACCUACUACCUCACCGGAUCCACAGUGUA

11. grep 提取gtf指定染色体注释信息

# 提取1号染色体gtf注释
zcat Homo_sapiens.GRCh38.112.gtf.gz|grep -w '^1'| \
gzip -c > Homo_sapiens.GRCh38.112.gtf.chr1.gz

12. GTF文件转换为bed格式文件

# 需要使用gtftobed, 需conda安装bedops
conda install bedops -y

# gtf转换为bed格式
# gzip -c 压缩为.gz格式
zcat Homo_sapiens.GRCh38.112.gtf.gz|gtf2bed --do-not-sort| \
gzip -c > Homo_sapiens.GRCh38.112.bed.gz

13. 提取GTF文件外显子exon信息

# 提取gtf exon信息
zcat Homo_sapiens.GRCh38.112.gtf.gz|grep exon | \
cut -f1,4,5,9 | cut -d ";" -f1  |  \
awk '{print $1, $2, $3, $5}' | \
sed -e 's/ /\t/g' | sed -e 's/\"//g' > gtf.exon.bed

微信扫一扫分享文章