24. 生信常见文件格式详解¶

>序列ID 描述信息（可选）
序列内容（可以折行，通常每行60或80个字符）

>NC_000001.11 Homo sapiens chromosome 1, GRCh38.p14
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC

>gi|556503834|ref|NC_000913.3| Escherichia coli K-12 MG1655
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTC
TGATAGCAGCTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGG

# ============ seqkit（推荐的FASTA/FASTQ瑞士军刀） ============

# 统计序列数量和总长度
seqkit stats genome.fasta
# 输出：file  format  type  num_seqs  sum_len  min_len  avg_len  max_len

# 提取指定ID的序列（-p 模式匹配）
seqkit grep -p "NC_000001" genome.fasta > chr1.fasta

# 获取前10条序列
seqkit head -n 10 genome.fasta > first10.fasta

# 将多行FASTA转为单行（方便grep处理）
seqkit seq -w 0 genome.fasta > oneline.fasta

# ============ samtools faidx（建索引+快速提取） ============

# 建FASTA索引（生成.fai文件）
samtools faidx genome.fasta

# 按区间提取序列（chr1的1000到2000位置）
samtools faidx genome.fasta chr1:1000-2000

@序列ID 描述信息
序列内容（一行，不折行）
+（分隔符，可选重复ID）
质量值字符串（和序列等长，每个字符对应一个碱基的质量）

@SRR1234567.1 1/1
NTGCACGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGATCGAT
+
#8ACCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
@SRR1234567.2 2/1
AGCTTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTA
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:

# ============ 基本统计 ============

# 快速统计（seqkit）
seqkit stats -a sample_R1.fastq.gz
# 输出包含：reads数、总碱基、Q20%、Q30%、GC含量

# 用zcat查看压缩FASTQ的前8行（2条reads）
zcat sample_R1.fastq.gz | head -8

# ============ 质控工具 ============

# FastQC 生成质量报告
fastqc -t 4 -o ./qc_report/ sample_R1.fastq.gz sample_R2.fastq.gz

# fastp 一步完成质控+过滤+报告（推荐）
fastp \
  -i raw_R1.fq.gz -I raw_R2.fq.gz \    # 输入：双端原始数据
  -o clean_R1.fq.gz -O clean_R2.fq.gz \ # 输出：过滤后数据
  -h report.html -j report.json \         # 报告文件
  --qualified_quality_phred 20 \          # Q20以下的碱基视为低质量
  --length_required 50 \                  # 过滤掉<50bp的reads
  --thread 4                              # 4线程

# ============ 进阶samtools操作 ============

# 查看BAM文件的header（包含参考基因组信息和程序记录）
samtools view -H sorted.bam

# 统计每条染色体的reads数
samtools idxstats sorted.bam
# 输出：染色体名  长度  mapped_reads  unmapped_reads

# 按区域提取（需要已排序+已建索引）
samtools view sorted.bam chr1:1000000-2000000 | head

# 计算测序深度（每个位置覆盖了多少reads）
samtools depth -a sorted.bam > depth.txt
# -a 表示输出所有位置（包括0深度）

# 统计比对质量分布
samtools view sorted.bam | awk '{print $5}' | sort | uniq -c | sort -rn

# ============ BAM转CRAM（节省存储） ============
samtools view -C -T reference.fasta -o output.cram sorted.bam
# -C 输出CRAM格式
# -T 指定参考基因组（CRAM必须）

# CRAM转BAM
samtools view -b -T reference.fasta -o output.bam input.cram

##fileformat=VCFv4.3
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1 Sample2
chr1    10177   rs367896724 A   AC  100 PASS    DP=500;AF=0.45  GT:DP:GQ    0/1:30:99   1/1:25:85
chr1    10235   .   T   TA  50  LowQual DP=100;AF=0.12  GT:DP:GQ    0/0:20:60   0/1:15:40
chr1    10352   rs201106462 T   TA  200 PASS    DP=800;AF=0.55  GT:DP:GQ    0/1:45:99   0/1:38:95

# ============ bcftools 常用操作 ============

# 查看VCF基本统计
bcftools stats input.vcf.gz > stats.txt

# 按质量过滤（QUAL>=30 且 深度>=10）
bcftools filter -i 'QUAL>=30 && INFO/DP>=10' input.vcf.gz -o filtered.vcf.gz

# GATK硬过滤推荐参数（SNP）
bcftools filter -i \
  'QD>=2.0 && FS<=60.0 && MQ>=40.0 && MQRankSum>=-12.5 && ReadPosRankSum>=-8.0' \
  snps.vcf.gz -o snps_filtered.vcf.gz

# 只提取PASS的变异
bcftools view -f PASS input.vcf.gz -o pass_only.vcf.gz

# 提取指定区域
bcftools view -r chr1:1000000-2000000 input.vcf.gz

# 提取指定样本
bcftools view -s Sample1,Sample2 input.vcf.gz -o subset.vcf.gz

# VCF转BCF（二进制压缩）
bcftools view -Ob input.vcf.gz -o output.bcf
bcftools index output.bcf

# 合并多个VCF
bcftools merge sample1.vcf.gz sample2.vcf.gz -o merged.vcf.gz

# 注释变异（添加rsID）
bcftools annotate -a dbsnp.vcf.gz -c ID input.vcf.gz -o annotated.vcf.gz

# 统计变异类型（SNP/InDel比例）
bcftools stats input.vcf.gz | grep "^SN"

##gff-version 3
chr1    ensembl gene    11869   14409   .   +   .   ID=ENSG00000223972;Name=DDX11L1;biotype=transcribed_unprocessed_pseudogene
chr1    ensembl mRNA    11869   14409   .   +   .   ID=ENST00000456328;Parent=ENSG00000223972;Name=DDX11L1-201
chr1    ensembl exon    11869   12227   .   +   .   ID=exon0001;Parent=ENST00000456328
chr1    ensembl exon    12613   12721   .   +   .   ID=exon0002;Parent=ENST00000456328
chr1    ensembl exon    13221   14409   .   +   .   ID=exon0003;Parent=ENST00000456328

chr1    ensembl gene    11869   14409   .   +   .   gene_id "ENSG00000223972"; gene_name "DDX11L1"; gene_biotype "transcribed_unprocessed_pseudogene";
chr1    ensembl transcript  11869   14409   .   +   .   gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; gene_name "DDX11L1";
chr1    ensembl exon    11869   12227   .   +   .   gene_id "ENSG00000223972"; transcript_id "ENST00000456328"; exon_number "1";

chr1    11868   14409   DDX11L1 0   +
chr1    14403   29570   WASH7P  0   -
chr1    29553   31109   MIR1302-2HG 0   +

# ============ 格式转换 ============

# GFF3转GTF（用gffread，来自Cufflinks/StringTie）
gffread annotation.gff3 -T -o annotation.gtf

# GTF转BED（用UCSC的gtfToGenePred + genePredToBed）
gtfToGenePred annotation.gtf /dev/stdout | genePredToBed /dev/stdin genes.bed

# 或用awk快速提取基因区间为BED
awk '$3=="gene"' annotation.gtf | \
  awk -F'\t' '{OFS="\t"; print $1, $4-1, $5, ".", "0", $7}' > genes.bed

# ============ bedtools 常用操作 ============

# 求两个BED文件的交集（找重叠区域）
bedtools intersect -a peaks.bed -b genes.bed > overlap.bed

# 求交集并报告重叠比例
bedtools intersect -a peaks.bed -b genes.bed -wo > overlap_with_size.bed

# 求不重叠的区域（A有但B没有）
bedtools subtract -a peaks.bed -b blacklist.bed > clean_peaks.bed

# 合并重叠或相邻的区间
bedtools merge -i sorted_peaks.bed > merged.bed

# 统计每个基因区域内有多少peaks
bedtools coverage -a genes.bed -b peaks.bed > gene_coverage.bed

# 获取BED区间对应的FASTA序列
bedtools getfasta -fi genome.fasta -bed regions.bed -fo regions.fa

# 在BED区间两侧各扩展1000bp
bedtools slop -i peaks.bed -g genome.sizes -b 1000 > expanded.bed

# 随机生成与输入BED相同大小的区间（用于统计背景）
bedtools shuffle -i peaks.bed -g genome.sizes > random_regions.bed

# 固定步长格式（fixedStep）
fixedStep chrom=chr1 start=1 step=10 span=10
10.5
12.3
8.7
15.2

# 可变步长格式（variableStep）
variableStep chrom=chr1 span=10
1000  15.3
1020  8.9
1050  22.1

# BAM转BigWig（用deeptools）
bamCoverage -b sorted.bam -o coverage.bw --normalizeUsing RPKM

# WIG转BigWig（需要chrom.sizes文件）
wigToBigWig input.wig chrom.sizes output.bw

# BigWig转BedGraph（可编辑的文本格式）
bigWigToBedGraph input.bw output.bedGraph

chr1    1000    2000    chr1    5000    6000    SV001   50  +   -
chr2    30000   31000   chr5    80000   81000   translocation1  99  +   +

Hugo_Symbol  Chromosome  Start_Position  End_Position  Variant_Classification  Variant_Type  Reference_Allele  Tumor_Seq_Allele2
TP53         chr17       7674220         7674220       Missense_Mutation       SNP           C                 T
BRCA1        chr17       43093488        43093488      Nonsense_Mutation       SNP           G                 A