import ibdutils.utils.ibdutils as ibdutils import ibdutils.runner.ibdne as ibdne # ---------------- prepare IBDNe input files ------------ # input files (simulated from single population model with selection) sp_ibd_files = [f"testdata/single_pop/tskibd/{i}.ibd" for i in range(1, 15)] sp_vcf_files = [f"testdata/single_pop/vcf/{i}.vcf.gz" for i in range(1, 15)] # parameters label_str = "tmp_sp_sim" ibdne_flatmeth = "none" ibdne_mincm = 2 # read ibd genome_14_100 = ibdutils.Genome.get_genome("simu_14chr_100cm") ibd = ibdutils.IBD(genome=genome_14_100, label=f"{label_str}_orig") ibd.read_ibd(ibd_fn_lst=sp_ibd_files) # remove highly related samples mat = ibd.make_ibd_matrix() unrelated_samples = ibd.get_unrelated_samples(ibdmat=mat) ibd.subset_ibd_by_samples(subset_samples=unrelated_samples) # prepare input for IBDNe # remove ibd with tmrca < 1.5 (required according to IBDNe paper) ibd.filter_ibd_by_time(min_tmrca=1.5) ibd.filter_ibd_by_length(min_seg_cm=ibdne_mincm) # calculate iHS ibd.calc_ihs(vcf_fn_lst=sp_vcf_files, min_maf=0.01) # calculate IBD coverage and find peaks ibd.calc_ibd_cov() ibd.find_peaks() # only keep peaks that contain a ihs hit ibd.filter_peaks_by_ihs() # plot IBD coverage with peaks marked with red shading ax = ibd.plot_coverage(which="ihsfilt") fig = ax.get_figure() fig.savefig("cov.png") # save IBD before remove peaks # of_orig_ibdne_obj = f"{label_str}_orig.ibdne.ibdobj.gz" # ibd.pickle_dump(of_orig_ibdne_obj) # saved ibd.obj file can be loaded by # ibd = ibdutils.IBD.pickle_load("xxx.ibd.obj.gz") # USEFUL information saved in the IBD object: ibd._df # IBD segment dataframe ibd._cov_df # IBD coverage dataframe ibd._peaks_df # IBD peaks and annotations # remove peaks ibd2 = ibd.duplicate(f"{label_str}_rmpeaks") ibd2.remove_peaks() ibd2._df = ibd2.cut_and_split_ibd() # USEFUL information saved in the IBD object: # IBD segment dataframe after peak removal ibd2._df # IBD coverage dataframe. Note this is not updated by `remove_peaks()` so now it # still reflects coverage before peak removal; if needed, you can call # `calc_ibd_cov()` to update _cov_df. ibd2._cov_df # IBD peaks and annotations. Note this is not updated by `remove_peaks()` so it # still reflects peaks before peak removal) ibd2._peaks_df # save IBD after remove peaks # of_rmpeaks_ibdne_obj = f"{label_str}_rmpeaks.ibdne.ibdobj.gz" # ibd2.pickle_dump(of_rmpeaks_ibdne_obj) ibdne_runner1 = ibdne.IbdNeRunner(ibd, ".", ".") ibdne_runner1.run(dry_run=True) # This generate three files *.ibd.gz/*.map/*.sh for running IBDNe ibdne_runner2 = ibdne.IbdNeRunner(ibd2, ".", ".") ibdne_runner2.run(dry_run=True) # This generate three files *.ibd.gz/*.map/*.sh for running IBDNe

import ibdutils.utils.ibdutils as ibdutils import ibdutils.runner.ibdne as ibdne import pandas as pd import numpy as np # ---------------- call infomap ------------ # input files mp_ibd_files = [f"testdata/multiple_pop/tskibd/{i}.ibd" for i in range(1, 15)] mp_vcf_files = [f"testdata/multiple_pop/vcf/{i}.vcf.gz" for i in range(1, 15)] # parameters label_str = "tmp_mp_sim" nsam = 200 # no. of isolates per subpopulation npop = 5 # no. of subpopulations transform = "square" # transformation of IBD matrix ntrials = 1000 # parameter of infomap algorithm # meta information with optional true labels (for purpose of comparison) meta = pd.DataFrame( { "Sample": np.arange(nsam * npop), # use haploid here "Population": np.repeat(np.arange(npop), nsam), } ) # read ibd from files genome_14_100 = ibdutils.Genome.get_genome("simu_14chr_100cm") ibd = ibdutils.IBD(genome=genome_14_100, label=f"{label_str}_orig") ibd.read_ibd(ibd_fn_lst=mp_ibd_files) # remove highly relatedness samples # NOTE: for empirical analysis, you might want to remove some short IBD segments # and mask ibd matrix elements of low values by changing arguments of the # `make_ibd_matrix` method, e.g.: # ibd.make_ibd_matrix(min_seg_cm=2, max_gw_ibd_cm=5.0) mat = ibd.make_ibd_matrix() unrelated_samples = ibd.get_unrelated_samples(ibdmat=mat) ibd.subset_ibd_by_samples(subset_samples=unrelated_samples) # calculate IBD coverage and find peaks ibd.calc_ibd_cov() ibd.find_peaks() # calculate iHS and filter peaks ibd.calc_ihs(vcf_fn_lst=mp_vcf_files, min_maf=0.01) ibd.filter_peaks_by_ihs() # plot IBD coverage with peaks marked with red shading ibd.plot_coverage(which="ihsfilt") # USEFUL information saved in the IBD object: ibd._df # IBD segment dataframe ibd._cov_df # IBD coverage dataframe ibd._peaks_df # IBD peaks and annotations # make a copy of the IBD object and remove IBD within peaks on the copy ibd2 = ibd.duplicate("rmpeak") ibd2.remove_peaks() # USEFUL information saved in the IBD object: # IBD segment dataframe after peak removal ibd2._df # IBD coverage dataframe. Note this is not updated by `remove_peaks()` so now it # still reflects coverage before peak removal; if needed, you can call # `calc_ibd_cov()` to update _cov_df. ibd2._cov_df # IBD peaks and annotations. Note this is not updated by `remove_peaks()` so it # still reflects peaks before peak removal) ibd2._peaks_df # run infomap on IBD object without peak removal mat = ibd.make_ibd_matrix() # run infomap on IBD object WITH peaks not removed member_df = ibd.call_infomap_get_member_df( mat, meta, trials=ntrials, transform=transform ) # run infomap on IBD object WITH peak removal mat2 = ibd2.make_ibd_matrix() member_df2 = ibd2.call_infomap_get_member_df( mat2, meta, trials=ntrials, transform=transform ) # Infomap identifies 2 main groups member_df.Rank.value_counts().iloc[:6] # 0 640 # 1 120 # 2 15 # 3 11 # 4 9 # 5 9 # Name: Rank, dtype: int64 # Infomap identifies 4 main groups when selection correction is performed member_df2.Rank.value_counts().iloc[:6] # 0 203 # 1 150 # 2 149 # 3 108 # 4 13 # 5 9 # Name: Rank, dtype: int64 # NOTE about Infomap result dataframe # column "Membership" contains the inferred labels (unsorted version) # column "Rank" contains the sorted version of inferred labels