云上管理基因型数据（以 DNANexus 平台、 UKB WES 为例）

1 相关软件

• KGGA 软件包（集成 GBC 、 CCF 基因组与大规模公共数据管理工具）

• Java 运行环境

  • 通过 wget https://download.oracle.com/java/24/latest/jdk-24_linux-x64_bin.tar.gz 下载 Oracle JDK

由于可能多次在不同的计算节点上使用 Java 环境，因此我们建议将相关资源存放在用户项目目录下。

2 实操流程

2.1 本地下载软件包

在本地设备上下载 KGGA，基因注释文件和 Oracle JDK：

# -c 启用断点续传
wget -c https://pmglab.top/kgga/download/lib/v1/kgga.jar
wget -c https://idc.biosino.org/pmglab/resource/kgg/kgga/resources/gene/refGene_hg38_kggseq_v2.txt.gz
wget -c https://download.oracle.com/java/24/latest/jdk-24_linux-x64_bin.tar.gz

2.2 上传到 RAP 平台

在 DNANexus 平台上进入项目路径，上传上述文件，获得文件的 ID （路径编号）：

KGGA: file-J21ZP4QJ5qQ25bjXZJX435Bq
JDK: file-J214ZzjJ5qQ8GfF272Fy4ZqP
refGene_hg38: file-J21Y8J0J5qQ6KZPxjvpgjF0P

2.3 获取数据路径

UKB 基因型数据储存在以下文件夹中：

# WES VCF
<PROJECT-ID>:/Bulk/Exome sequences/Population level exome OQFE variants, pVCF format - final release/

# WGS VCF
<PROJECT-ID>:/Bulk/DRAGEN WGS/ML-corrected DRAGEN population level WGS variants, pVCF format [500k release]/

# WGS PGEN
<PROJECT-ID>:/Bulk/DRAGEN WGS/DRAGEN population level WGS variants, PLINK format [500k release]/

有关基因型的存储说明，请参考：https://biobank.ndph.ox.ac.uk/ukb/label.cgi?id=100314

所有文件都以 ukb23157_c<CHROM>_b<BLOCK>_v1.vcf.gz 命名。

2.4 开启计算节点并配置运行环境

计算节点配置: mem1_ssd1_v2_x4，7.8 GB 内存，93 GB 硬盘，4 核心。

创建 Terminal 工作区:

cd /opt/notebooks/
dx download file-J21ZP4QJ5qQ25bjXZJX435Bq
dx download file-J214ZzjJ5qQ8GfF272Fy4ZqP
tar zxvf ./jdk-24_linux-x64_bin.tar.gz
export PATH=/opt/notebooks/jdk-24.0.2/bin:$PATH

2.5 下载数据到计算节点

由于 UKB 单个 VCF 文件很大，无法一次性加载所有的子文件进行处理。因此，我们先将整个染色体的数据均等分为 $k$ 份，然后在每一个节点上处理其中的一份数据（一份数据包含多个此染色体的子文件）。处理过程包括：基因型水平、个体水平质控，格式转换。处理完成的 GTB 文件回传到 DNANexus 平台。最后，整个染色体的数据处理完成时，再到计算节点上快速拼接为完整的 GTB 文件。

# ---------- 任务基本信息 ----------
  CHROM="21"
  # 总批次数（将对应染色体下的文件分为多少份）
  BATCH_NUM=5  
  # 当前批次索引（从0开始，范围 0~BATCH_NUM-1）
  BATCH_INDEX=0
  # 输入数据路径
  INPUT_PROJECT_ID="J04yP4QJ5qQ4pgP2xQpb0Qp3"
  INPUT_PROJECT_DIR="/Bulk/Exome sequences/Population level exome OQFE variants, pVCF format - final release/"
  # 输出数据路径, 【请确保此文件夹存在，否则会在 upload 步骤出错】
  OUTPUT_PROJECT_ID="J04yP4QJ5qQ4pgP2xQpb0Qp3"
  OUTPUT_PROJECT_DIR="/user/zlb/ukb/wes/chr${CHROM}"

# ---------- 文件列表 ----------
  echo "[STEP 1] 查找染色体 ${CHROM} 的 VCF 文件..."
  # 搜索 ${INPUT_PROJECT_ID}:${INPUT_PROJECT_DIR} 下的所有 vcf.gz 文件
  mapfile -t FILES < <(
    dx find data -r --name "ukb23157_c${CHROM}_b*.vcf.gz" \
      --path "project-${INPUT_PROJECT_ID}:${INPUT_PROJECT_DIR}/" \
      --brief \
    | xargs -n1 -I {} dx describe {} --name
  )

  TOTAL_FILES=${#FILES[@]}
  echo "发现 ${TOTAL_FILES} 个文件"

# ---------- 任务映射 ----------
  if [[ $BATCH_NUM -le 0 ]]; then
    echo "错误：BATCH_NUM 必须大于0"
    exit 1
  fi

  # 计算每个批次的大小（向上取整）
  BATCH_SIZE=$(( (TOTAL_FILES + BATCH_NUM - 1) / BATCH_NUM ))
  START_INDEX=$(( BATCH_INDEX * BATCH_SIZE ))
  END_INDEX=$(( START_INDEX + BATCH_SIZE - 1 ))

  # 校验任务批次索引合法性
  if [[ $BATCH_INDEX -ge $BATCH_NUM ]]; then
    echo "错误：BATCH_INDEX 不能超过 BATCH_NUM-1"
    exit 1
  fi

  # 处理最后一个批次的越界情况
  if [[ $END_INDEX -ge $TOTAL_FILES ]]; then
    END_INDEX=$(( TOTAL_FILES - 1 ))
  fi

  # 空批次检查
  if [[ $START_INDEX -gt $END_INDEX ]]; then
    echo "警告：当前批次为空（文件已处理完毕）"
    exit 0
  fi

  echo "总批次数：${BATCH_NUM} ｜ 每批约处理 ${BATCH_SIZE} 个文件"
  echo "当前批次：${BATCH_INDEX} (文件范围 ${START_INDEX}-${END_INDEX})"

# ---------- 处理程序 ----------
  for (( i=START_INDEX; i<=END_INDEX; i++ ))
  do
      INPUT="${FILES[$i]}"
      echo "进度：$((i - START_INDEX + 1))/$((END_INDEX - START_INDEX + 1)) - 文件 ${INPUT}"

      # 输出文件, 将 ukb23157_c1_b1.vcf.gz 改名为 ukb23157_c1_b1.gtb
      OUTPUT="$(basename ${INPUT} .vcf.gz).gtb"

      # 检查输出文件是否已经存在，如果存在则跳过
      if dx ls "project-${OUTPUT_PROJECT_ID}:${OUTPUT_PROJECT_DIR}/${OUTPUT}" > /dev/null 2>&1; then
          echo "文件已存在，跳过: ${OUTPUT}"
          continue
      fi

      # 下载文件
      if ls "${INPUT}" > /dev/null 2>&1; then
          echo "本地文件已存在，跳过下载: ${INPUT}"
      else
          echo "下载文件: ${INPUT}"
          dx download "project-${INPUT_PROJECT_ID}:${INPUT_PROJECT_DIR}/${INPUT}"
      fi

      # 处理文件
      echo "转换 VCF -> GTB..."
      java -jar -Xms6g -Xmx6g kgga.jar gbc convert \
           vcf2gtb "${INPUT}" -o "${OUTPUT}" \
           --field --prune --seq-an 1~ --seq-af 0.000001~0.999999 \
           -t 4 || { echo "Java conversion failed for ${INPUT}"; exit 1; }

      # 上传文件
      echo "上传文件: ${OUTPUT}"
      dx upload "${OUTPUT}" --path "project-${OUTPUT_PROJECT_ID}:${OUTPUT_PROJECT_DIR}/${OUTPUT}" || { echo "Upload failed for ${OUTPUT}"; exit 1; }

      # 删除处理好的本地文件
      rm "${INPUT}" "${OUTPUT}" || { echo "Failed to clean up files"; exit 1; }
      echo "处理完成: ${INPUT}"
  done

echo "批次 ${BATCH_INDEX} 处理完成"

同一个染色体的所有数据处理完成后，执行以下指令合并文件：

# ---------- 任务基本信息 ----------
  CHROM="21"
  # 输入数据路径
  INPUT_PROJECT_ID="J04yP4QJ5qQ4pgP2xQpb0Qp3"
  INPUT_PROJECT_DIR="/.table-exporter/user/UKB-WES GTB Format/chr${CHROM}"
  # 输出数据路径
  OUTPUT_PROJECT_ID="J04yP4QJ5qQ4pgP2xQpb0Qp3"
  OUTPUT_PROJECT_DIR="/user/zlb/ukb/wes/"
  OUTPUT="chr${CHROM}.hg38.gtb"

# ---------- 处理程序 ----------
  # 获取所有文件列表
  mapfile -t FILES < <(
    dx find data --name "ukb23157_c${CHROM}_b*.gtb" \
      --path "project-${INPUT_PROJECT_ID}:${INPUT_PROJECT_DIR}/" \
      --brief \
    | xargs -n1 -I {} dx describe {} --name
  )

  for (( i=0; i<=${#FILES[@]}; i++ ))
  do
      dx download "project-${INPUT_PROJECT_ID}:${INPUT_PROJECT_DIR}/${FILES[$i]}"
  done

  # 按照编号顺序排序
  SORTED_FILES=($(ls ./*.gtb | sort -t '_' -k3.2n))

  # 拼接 GTB 文件
  java -jar kgga.jar gbc concat "${SORTED_FILES[@]}" -o "${OUTPUT}"

  # 回传数据
  dx upload "${OUTPUT}" --path "project-${OUTPUT_PROJECT_ID}:${OUTPUT_PROJECT_DIR}/${OUTPUT}" || { echo "Upload failed for ${OUTPUT}"; exit 1; }

  # 删除子文件
  rm -rf ukb23157_c*.gtb

由此形成单个单个染色体的 GTB 文件。

2.6 从 GTB 转换到特定的格式 (基因型水平分析)

GTB 格式相比其他存储格式通常能够节省 50%~99% 的存储开销，而其分析和计算时的性能能够媲美高性能的位编码计算格式（PLINK 的PGEN 和 BED 格式），简易文档见 https://pmglab.top/kgga/v1.0/GBCQuickStart.html。

访问特定的子样本和子区间数据，并导出为 PGEN 格式 （先使用 pip3 install jep zstandard pgenlib 配置依赖）：

java -Djava.library.path="$(pip3 show jep | grep Location | awk "{print \$2\"/jep\"}")" -jar \
    kgga.jar gbc convert gtb2plink chr21.hg38.gtb --seq-af 0.05~0.95 -ot pgen

2.7 基于GTB用KGGA对位点进行基因属性的注释，只保留影响蛋白质的非同义突变

基于 GTB 格式， KGGA 可以对位点进行高速基因组属性注释，本范例只演示对位点进行基因属性的注释和过滤。 更详细的功能见 KGGA 的在线文档。

位点进行基因属性的注释和过滤，注释结果保存到 TSV文本文件中：

#准备KGGA注释用的资源路径
mkdir resources
cd resources
mkdir gene
cd gene

#下载KGGA注释用的资源文件
dx download file-J21Y8J0J5qQ6KZPxjvpgjF0P
cd ../../

#正式开始注释
   java -Xmx4g -jar kgga.jar annotate \
   --input-gty-file ./chr21.hg38.gtb \
   --output ./annot \
   --seq-ac 1 \
   --gene-model-database refgene \
   --gene-feature-included 0~6

通过Linux 命令快速查看注释结果

zless /opt/notebooks/annot/variants.hg38.tsv.gz