diff --git a/build_pop.py b/build_pop.py
index 356fe9837c33893f19fa77a4af9ac246c4b55ce7..eae3f35ccddb624a849f900242a2c8902b4549e8 100755
--- a/build_pop.py
+++ b/build_pop.py
@@ -38,6 +38,7 @@ def parse_args():
                         type=int, default=30)
     parser.add_argument("-q", "--quiet", help="Don't ask anything, choose default answer instead", action="store_const",
                         const=True, default=False)
+    parser.add_argument("-t", "--threads", help="Number of threads", default=4)
 
     args = parser.parse_args()
     return args
@@ -260,7 +261,7 @@ def build_fastas_chromosomes(reference, genotypes_for_inds, haploid, output_dir)
 
 
 def generate_samples_fastq(haploid, nb_inds, output_dir, coverage, read_len, insert_len_mean, insert_len_sd,
-                           prg_path):
+                           prg_path, threads):
 
     print("GENERATE RANDOM READS FOR EACH INDIVIDUAL FROM GENOME...\n")
 
@@ -268,21 +269,21 @@ def generate_samples_fastq(haploid, nb_inds, output_dir, coverage, read_len, ins
     if not haploid:
         cmd = str("{4}pirs/pirs simulate -z -x {3} -d -B {4}pirs/Profiles/Base-Calling_Profiles/humNew.PE100.matrix.gz "
                   "-I {4}pirs/Profiles/InDel_Profiles/phixv2.InDel.matrix -l {5} -m {6} -v {7} "
-                  "-G {4}pirs/Profiles/GC-depth_Profiles/humNew.gcdep_100.dat -o {0} {1} {2}")
+                  "-G {4}pirs/Profiles/GC-depth_Profiles/humNew.gcdep_100.dat -o {0} -t {8} {1} {2}")
     else:
         cmd = str("{4}pirs/pirs simulate -z -x {3} -B {4}pirs/Profiles/Base-Calling_Profiles/humNew.PE100.matrix.gz "
                   "-I {4}pirs/Profiles/InDel_Profiles/phixv2.InDel.matrix -l {5} -m {6} -v {7} "
-                  "-G {4}pirs/Profiles/GC-depth_Profiles/humNew.gcdep_100.dat -o {0} {1}")
+                  "-G {4}pirs/Profiles/GC-depth_Profiles/humNew.gcdep_100.dat -o {0} -t {8} {1}")
     for i in range(1, nb_inds+1):
         prefix = os.path.join(output_dir, "INDIV_" + str(i))
         chr0 = prefix + "_chr_0.fasta"
         chr1 = prefix + "_chr_1.fasta"
         if not haploid:
             os.system(cmd.format(prefix, chr0, chr1, coverage, prg_path + os.path.sep, read_len,
-                                 insert_len_mean, insert_len_sd))
+                                 insert_len_mean, insert_len_sd, threads))
         else:
             os.system(cmd.format(prefix, chr0, "", coverage, prg_path + os.path.sep, read_len,
-                                 insert_len_mean, insert_len_sd))
+                                 insert_len_mean, insert_len_sd, threads))
 
 
 def confirm(deletions: dict):
@@ -336,7 +337,7 @@ def main():
                                                       tmp_dir, prg_path)
         build_fastas_chromosomes(reference, genotypes_for_inds, haploid, output_dir)
         generate_samples_fastq(haploid, nb_inds, output_dir, args.coverage, args.read_len, args.insert_len_mean,
-                               args.insert_len_sd, prg_path)
+                               args.insert_len_sd, prg_path, args.threads)
         print("DONE!\n")
     else:
         print("Aborted!\n")