Unconfigured Ad

Collapse
X
 
  • Time
  • Show
Clear All
new posts
  • Lien
    Member
    • Dec 2009
    • 47

    #1

    Question about python script generating synthetic sequences

    Hi all,

    I'm trying to replicate a Nature Protocols paper: genome-wide copy number analysis of single cells from Baslan et al., 2012.
    In step 69 they generate synthetic 50 bp reads for each position on the hg19 genome: start at position one in the chromosome sequence and take the first 50 bases. Also create read ID strings and quality score strings. These can be Illumina or fastq format. Output these to a file and continue likewise at positions 2, 3, 4 and so on until the end of the chromosome is reached.
    The script to do this can be found in the supplementary data, and I've pasted it below:

    #!/usr/bin/env python

    import re
    import sys
    import random
    import string

    def main():

    CHROMLIST = open("/filepath/chromlist.txt", "r")

    chromlist = []
    for x in CHROMLIST:
    if x.rstrip().find("_") == -1:
    chromlist.append(x.rstrip())

    INFILE = False
    OUTFILE = open("/filepath/sequence.part.0.k50.txt", "w")
    outfilecount = 0
    counter = 0

    for chrFile in chromlist:
    a = chrFile.split("/")
    b = a[len(a) - 1]
    thisChrom = b.split(".")[0]

    if INFILE:
    INFILE.close()
    INFILE = open(chrFile, "r")
    chr = []
    x = ""
    y = ""
    INFILE.readline()
    for x in INFILE:
    chr.append(x.rstrip())
    x = "".join(chr)
    y = x.upper()
    print "after read " + thisChrom
    sys.stdout.flush()

    for i in range(len(y) - 50):
    counter += 1
    thisRead = y[i: (i+50)]

    if counter % 150000000 == 0:
    OUTFILE.close()
    outfilecount += 1
    OUTFILE = open("/filepath/sequence.part." + str(outfilecount) + ".k50.txt", "w")

    OUTFILE.write("@" + thisChrom + "." + str(i))
    OUTFILE.write("\n")
    OUTFILE.write(thisRead)
    OUTFILE.write("\n")
    OUTFILE.write("+" + thisChrom + "." + str(i))
    OUTFILE.write("\n")
    OUTFILE.write("b" * len(thisRead))
    OUTFILE.write("\n")

    INFILE.close()

    OUTFILE.close()


    if __name__ == "__main__":
    main()


    I don't have a lot of experience with python, and when I run this script (after adapting the filepath), I only get 1 empty file. Is there something else I need to adapt?

    Many thanks in advance,
    Lien
    Last edited by GenoMax; 04-18-2016, 02:57 AM. Reason: Moved post to a new thread to increase visibility
  • GenoMax
    Senior Member
    • Feb 2008
    • 7142

    #2
    For reference cross-posted: https://www.biostars.org/p/187267/

    Comment

    Latest Articles

    Collapse

    • SEQadmin2
      Multiomic Liquid Biopsy: What DNA/RNA Panels Reveal That DNA Alone Can’t
      by SEQadmin2



      Tumors contain complex mixtures of cell populations that evolve in response to treatment and host-specific factors. Tissue biopsies provide only static snapshots of a tumor’s state. But what if we could follow that evolution over time using a minimally invasive technique?

      Liquid biopsies can detect tumor-derived material circulating in the bloodstream, providing a way to study treatment response, detect emerging resistance, and discover new biomarkers. Circulating
      ...
      10-14-2026, 09:55 AM
    • SEQadmin2
      New Genomics Technologies Take Aim at Long-Standing Limits
      by SEQadmin2


      Researchers using sequencing and genomics tools often have to make trade-offs. They can choose between speed or scale, short reads or long-range information, or targeted panels or a view of the whole transcriptome. New technologies that have been released this year are built to address those tough choices.

      We asked six companies the same four questions to learn about their latest products. The new technologies bring a lot to the table, including rethinking sequencing
      ...
      09-28-2026, 10:25 AM

    ad_right_rmr

    Collapse

    News

    Collapse

    Topics Statistics Last Post
    Started by SEQadmin2, Yesterday, 11:13 AM
    0 responses
    15 views
    0 reactions
    Last Post SEQadmin2  
    Started by SEQadmin2, Yesterday, 10:48 AM
    0 responses
    15 views
    0 reactions
    Last Post SEQadmin2  
    Started by SEQadmin2, 10-15-2026, 10:52 AM
    0 responses
    35 views
    0 reactions
    Last Post SEQadmin2  
    Started by SEQadmin2, 10-14-2026, 05:39 AM
    0 responses
    42 views
    0 reactions
    Last Post SEQadmin2  
    Working...