生物信息文件夹

NIPT麻木再来一遍

  • 2022年 10月21日

前文再续，书接上回。

重新梳理一遍NIPT的分析方法。

QC

直接使用fastp对原始数据进行质控，原始数据是SE50。

fastp -i sample.fastq.gz -o sample.clean.fq.gz -w 16 -j sample.json -h sample.html

比对

这一步有用bwa mem、bwa aln和bowtie2的，差异不是很大。但是，使用的参考基因组会有影响，例如使用hs37d5比ucsc.hg19稍好，因为hs37d5会预先mask掉PAR区域。比对完后要去重。

bwa mem -t 32 hs37d5.fa sample.clean.fq.gz | samtools view -bSh - | samtools sort -@ 32 - -o sample.sort.bam
gatk MarkDuplicates -I sample.sort.bam -O sample.rmdup.bam -M sample.dups.txt --REMOVE_DUPLICATES true

Bed

需要对参考基因组进行划窗，用于GC校正。

PAR区

Version Region Note
GRCh38 chrY:10001-2781479/chrX:10001-2781479 PAR1
GRCh38 chrY:56887903-57217415/chrX:155701383-156030895 PAR2
GRCh37 chrY:10001-2649520/chrX:60001-2699520 PAR1
GRCh37 chrY:59034050-59363566/chrX:154931044-155260560 PAR2

使用bedtools，根据广东省的团标，划窗可以是20kb、60kb、100kb。这里使用50kb。

可以使用fai索引文件来划分

bedtools makewindows -g hs37d5.fa.fai -w 50000 > hs37d5.50kb.bed

也可以使用预先建立好的bed进行划分

bedtools makewindows -b hs37d5.bed -w 50000 > hs37d5.50kb.bed

计算各个窗口的GC比例

bedtools nuc -fi hs37d5.fa -bed hs37d5.50kb.bed > hs37d5.50kb.gc.bed

去掉与PAR区有重合区域的窗口（可选）

bedtools intersect -a hs37d5.50kb.gc.bed -b PAR.bed -v > hs37d5.50kb.gc.noPAR.bed

去掉GC比例过低和过高的区域（可选），注意自己的bed中GC在第几列。

cat hs37d5.50kb.gc.noPAR.bed | awk '{if(($6>=0.2 && $6<=0.8)) print $0}' > hs37d5.50kb.gc.noPAR.filter.bed

去掉黑名单区域（可选）

bedtools intersect -a hs37d5.50kb.gc.noPAR.filter.bed -b wgEncodeDukeMapabilityRegionsExcludable.bed > hs37d5.50kb.gc.noPAR.filter.rmBlack.bed

统计

统计各个窗口的read counts，这会将counts补充到最后一列

bedtools intersect -a hs37d5.50kb.gc.noPAR.filter.rmBlack.bed -b sample.rmdup.bam -c -wa > sample.counts.txt

GC校正

使用loess回归校正read counts，代码是R

args <- commandArgs(T)
rawCounts <- args[1]
corCounts <- args[2]

# V13是counts，V5是GC，按实际修改
RC_DT <- read.table(rawCounts, sep="\t", head=FALSE)
gcCount.loess <- loess(
  V13~V5,
  data=RC_DT,
  control=loess.control(surface="direct"),
  degree=2
)
predictions <- predict(gcCount.loess, RC_DT$V5)

RC_DT$RC <- predictions
RC_DT$RCor <- ifelse(RC_DT$RC<0, 0, round(RC_DT$RC, 0))

write.table(RC_DT, corCounts, sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE)

胎儿浓度

男胎浓度

接下来需要计算胎儿浓度。易知，男胎浓度使用chrY浓度来代替即可，但由于女胎也有少量reads会比对到chrY上，因此需要一定的公式校正。这里参考defrag算法。

代码为python

import pandas as pd

# 把校正后UR代入，再使用UR均值校正一次UR
# URAdj = UR * (URmean / URGC)

def GetURDateFrame(corFile):
    df = pd.read_csv(corFile, header=None, sep="\t")
    # df.rename(columns={0: "chrom", 1: "start", 2: "end", 3: "name", 5: "GC", 13: "UR", 15: "URGC"}, inplace=True)
    df.rename(columns={0: "chrom", 1: "start", 2: "end", 4: "GC", 12: "UR", 14: "URGC"}, inplace=True)
    # df2 = df[["chrom", "start", "end", "name", "GC", "UR", "URGC"]]
    df2 = df[["chrom", "start", "end", "GC", "UR", "URGC"]]

    # UR均值
    URmean = df2["UR"].mean()

    # UR修正值
    df2.insert(loc=6, column="URAdj", value=0)
    nonZeroList = df2.loc[df2["URGC"] != 0].index
    for i in nonZeroList:
        URGC = df2.loc[i, "URGC"]
        UR = df2.loc[i, "UR"]
        URAdj = round(UR * (URmean / URGC))
        # 还有一种校正方式
        # URAdj = round(UR - (URGC - URmean))
        df2.loc[i, "URAdj"] = URAdj
    return df2

defrag算法计算男胎胎儿浓度，即用（待测样本的Y丰度 - 若干女胎样本的平均Y丰度）除以（若干男性成人平均Y丰度 - 若干女胎样本的平均Y丰度）

\[FF = \frac{\%Y_{XY}-\%Y_{XX}}{\%Y_{man}-\%Y_{XX}}\]

在不改变实验环境与条件的前提下，$\%Y_{man}$和$\%Y_{XX}$是两个常数。

# 获得Y丰度
def percY(dataFrame):
    sumY = dataFrame[dataFrame["chrom"] == "chrY"]["URAdj"].sum()
    sumAll = dataFrame["URAdj"].sum()
    Y_per = float(sumY) / sumAll
    return Y_per

# 获得一批次的平均Y丰度
def meanPerY(dir):
    Y_per_List = []
    for i in os.listdir(dir):
        df = GetURDateFrame(dir + "/" + i)
        y_p = percY(df)
        Y_per_List.append(y_p)
    meanYPer = sum(Y_per_List) / len(Y_per_List)
    return meanYPer

女胎浓度

和男胎不同，女胎浓度不能通过chrX的丰度获得。seqFF算法根据常染色体read counts（去掉13，18，21三个染色体）和男胎浓度来训练获得模型，再把模型用于预测胎儿浓度，模型同时适用于女胎。模型是同时使用弹性网络模型获得预测值Ye，以及分级加权模型获得预测值Yw，最终预测胎儿浓度为Ye和Yw的平均值。seqFF的脚本可以在这里下载到，已经预先训练了模型。想尝试一下怎么把模型训练出来和以及怎么获得那两个校正常数的，但是脚本中没有相关信息。

将counts处理为如下格式的矩阵

checkPoint sample1 sample2 sample3
chr1-50000 74 97 188
chr1-100000 57 83 207
chr1-150000 63 63 141

弹性网络的sklearn

import pandas as pd
from sklearn.linear_model import ElasticNetCV
import joblib

# 训练
df = pd.read_csv("train_male.txt", sep="\t", header=0, index_col=0)
dft = df.T
cols = [i for i in dft.columns if i != "ff"]
X = dft[cols]
y = dft["ff"]
regr = ElasticNetCV(cv=10, random_state=0, max_iter=1000)
regr.fit(X, y)
## 输出模型，可以复用
joblib.dump(regr, "output.model")

# 测试
## 载入模型
regr = joblib.load("output.model")
df_test = pd.read_csv("test_female.txt", sep="\t", header=0, index_col=0)
dft_test = df_test.T
cols_test = [i for i in dft_test.columns if i != "ff"]
X_test = dft_test[cols_test]
y_test_act = list(dft_test["ff"])
predict_test = list(regr.predict(X_test))
zipped = list(zip(y_test_act, predict_test))
for z in zipped:
    print(z[0], z[1])

WRSC没看懂怎么做，摸索中。

Z值计算

Z值计算算法比较简单，参考过往此篇。

搜索