fastq文件单个和批量转换至bigwig文件

一般我们可以直接从NCBI上下载基因的fastq文件。

​​FASTQ文件​​是存储生物测序原始数据的标准文本格式,主要用于记录测序得到的​​核酸序列(如DNA/RNA)及其对应的测序质量评分​​。它是生物信息分析(如基因组组装、变异检测等)的基础数据格式。
1.数据比对:FASTQ → BAM

bowtie2-build reference.fasta reference_index  构建参考基因组索引(仅需一次)

[1] bowtie2 -x reference_index -1 SRR8268703_1.fastq -2 SRR8268703_2.fastq -S output.sam --threads 8  # 双端测序数据比对

[2] bowtie2 -x reference_index -U SRR8268730.fastq -S output.sam --threads 8  # 单端测序数据比对

samtools sort -@ 8 -o output.bam output.sam
samtools index output.bam  # 建立索引

2.生成覆盖度文件:BAM → BedGraph

bedtools genomecov -ibam output.bam -bga -split > output.coverage.bedgraph  # 使用bedtools生成覆盖度文件

3.转换BedGraph为BigWig

bedGraphToBigWig output.coverage.bedgraph chrom.sizes output.bw 

4.从NCBI上批量下载的基因一般在一个文件夹下,可以使用fastq_to_bw.sh脚本进行批量转换。

首先找到项目数据的登录号在“SRA”数据库中搜索。随便选中一个SRR号,然后点击“All runs”。接着选中你想下载的数据,点击"Accession list",会下载一个包含选中数据SRR号的文件(SRR_Acc_List.txt)。

接着可以进行批量下载:

prefetch -O output --option-file SRR_Acc_List.txt

我一般会保存在./sra-data里,如图:

将sra文件转化为fastq格式:

ls .sra | parallel -j 4 "fasterq-dump --split-files --outdir ./fastq {} && gzip ./fastq/{}.fastq"

但是大部分时间我们想直接获取fastq文件,则可以:

cat SRR_Acc_List.txt | parallel -j 4 "fasterq-dump --progress --split-files --outdir ./fastq {}"  # 将fastq文件保存至当前文件夹的fastq文件夹中

之后就可以用fastq_to_bw.sh脚本了~

#!/bin/bash
set -euo pipefail # 出现错误时退出

# 配置参数
REF_INDEX="reference_index" # 参考基因组索引前缀
CHROM_SIZES="chrom.sizes"   # 染色体尺寸文件
FASTQ_DIR="../fastq"        # FASTQ文件所在目录
THREADS=${LSB_DJOB_NUMPROC:-4} # 默认使用4个线程或LSF调度系统分配的线程数

# 创建输出目录
mkdir -p aligned_bam bedgraph bigwig logs

# 自动识别所有样本(支持双端和单端)
find "$FASTQ_DIR" -maxdepth 1 -name "*.fastq" \
  | sed -E 's/(_1|_2)?\.fastq$//' \
  | sort -u \
  | while read SAMPLE_PREFIX; do
    SAMPLE=$(basename "$SAMPLE_PREFIX")
    R1="${FASTQ_DIR}/${SAMPLE}_1.fastq"
    R2="${FASTQ_DIR}/${SAMPLE}_2.fastq"
    SINGLE="${FASTQ_DIR}/${SAMPLE}.fastq"

echo "====== 处理样本 $SAMPLE ======" | tee -a logs/${SAMPLE}.log

if [[ -f "$R1" && -f "$R2" ]]; then
    echo "[INFO] 双端测序" | tee -a logs/${SAMPLE}.log
    CMD="bowtie2 -x $REF_INDEX -1 $R1 -2 $R2 --threads $THREADS"
elif [[ -f "$SINGLE" ]]; then
    echo "[INFO] 单端测序" | tee -a logs/${SAMPLE}.log
    CMD="bowtie2 -x $REF_INDEX -U $SINGLE --threads $THREADS"
else
    echo "[ERROR] 文件不存在: $SAMPLE" | tee -a logs/${SAMPLE}.log
    exit 1
fi

# 比对和排序
echo "[INFO] 开始比对..." | tee -a logs/${SAMPLE}.log
$CMD \\
    | samtools view -@ $THREADS -bS - \\
    | samtools sort -@ $THREADS -o aligned_bam/${SAMPLE}_sorted.bam 2>> logs/${SAMPLE}.log

# 生成索引
echo "[INFO] 生成BAM索引..." | tee -a logs/${SAMPLE}.log
samtools index aligned_bam/${SAMPLE}_sorted.bam 2>> logs/${SAMPLE}.log

# 生成BedGraph
echo "[INFO] 生成覆盖度文件..." | tee -a logs/${SAMPLE}.log
bedtools genomecov -ibam aligned_bam/${SAMPLE}_sorted.bam -bg > bedgraph/${SAMPLE}.bedgraph 2>> logs/${SAMPLE}.log

# 转换为BigWig
echo "[INFO] 转换到BigWig格式..." | tee -a logs/${SAMPLE}.log
bedGraphToBigWig bedgraph/${SAMPLE}.bedgraph $CHROM_SIZES bigwig/${SAMPLE}.bw 2>> logs/${SAMPLE}.log

echo "[INFO] 处理完成!生成文件:bigwig/${SAMPLE}.bw" | tee -a logs/${SAMPLE}.log
done

echo "====== 全部样本处理完成 ======"

评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值