一般我们可以直接从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 "====== 全部样本处理完成 ======"