Untitled

mail@pastecode.io avatar
unknown
plain_text
a year ago
1.1 kB
5
Indexable
### Filtering variants ###
# Reading HapMap3+ variants
info <- readRDS("/home/chronos/trp/map_hm3_plus.rds")
    
# Reading pathway info
bim_path <- "/home/chronos/tdora_uj/sirt1_data/new/comorbid_pathway.bim"
pathway_info <- fread(bim_path, data.table = FALSE)  # reading bim file
pathway_info <- setNames(pathway_info[c(1, 4)], c("chr", "pos"))

# Merging the pathway to HapMap3+ variants
info_plus <- full_join(
    info[c("chr", "pos")], pathway_info, by = c("chr", "pos")
)

# Filtering by position
filtered_sumstats <- inner_join(sumstats, info_plus, by = c("chr", "pos"))

### Matching discovery and target genetic data ###
info_snp <- snp_match(
    filtered_sumstats, map, join_by_pos = TRUE, return_flip_and_rev = TRUE
)

# Creating a subset for the pathway of interest
pathway_snp <- inner_join(info_snp, pathway_info, by = c("chr", "pos"))
pathway_indices <- which(df_beta[["_NUM_ID_"]] %in% pathway_snp[["_NUM_ID_"]])

### Run LDpred2-auto as usual ###

### Pathway-based polygenic risk ###
pathway_auto <- big_prodVec(
    G, beta_auto[pathway_indices], ind.col = pathway_snp[["_NUM_ID_"]]
)