Cart

×

You have 0 items in your cart.

# Code for chapter 18

#### For instructors:

Below are listed all the scripts as shown in chapter 18 of the book. These are free to copy and paste into your code editor.

## Script 18.1

 script_18_1.py # Script 18.1mass = 3.4volume = 1.8density = mass/volume       # Divisiondensity = round(density, 3) # Round to three decimal placesprint(density) script_18_1.out [ah@hobbes Documents]\$ python script_18_1.py 1.889[ah@hobbes Documents]\$

## Script 18.2

 script_18_2.py # Script 18.2x = 2     # Integery = 5     # Integerz = 3.0   # Floating pointprint(x * y) # 10  - Integerprint(x / y) # 0.4 - Floating point due to divisionprint(x + z) # 5.0 - Floating pointt = type(x * y - z)  # Get a value's data typeprint(t)             # float script_18_2.out [ah@hobbes Documents]\$ python script_18_2.py 1005.0[ah@hobbes Documents]\$

## Script 18.3

 script_18_3.py # Script 18.3x = 1            # An arbitrary integerprint(x.__doc__) # Print Python documentation for this object script_18_3.out [ah@hobbes Documents]\$ python script_18_3.py int(x=0) -> int or longint(x, base=10) -> int or longConvert a number or string to an integer, or return 0 if no argumentsare given.  If x is floating point, the conversion truncates towards zero.If x is outside the integer range, the function returns a long instead.If x is not a number or if base is given, then x must be a string orUnicode object representing an integer literal in the given base.  Theliteral can be preceded by '+' or '-' and be surrounded by whitespace.The base defaults to 10.  Valid bases are 0 and 2-36.  Base 0 means tointerpret the base from the string as an integer literal.>>> int('0b100', base=0)4[ah@hobbes Documents]\$

## Script 18.4

 script_18_4.py # Script 18.4x = Trueprint(x)print(bool(0.0)) # Falseprint(bool(7))   # True script_18_4.out [ah@hobbes Documents]\$ python script_18_4.py TrueFalseTrue[ah@hobbes Documents]\$

## Script 18.5

 script_18_5.py # Script 18.5x = 2print(x > 5) # Falsey = x < 3print(y)     # True script_18_5.out [ah@hobbes Documents]\$ python script_18_5.py FalseTrue[ah@hobbes Documents]\$

## Script 18.6

 script_18_6.py # Script 18.6x = 1y = -1x > y and y > 1   # False - second comparison is Falseprint(x > y and y > 1)x != 2 or y == 1  # True - first comparison is True; x does not equal 2print(x != 2 or y == 1)not (y > x)       # True - comparison is False; not False is Trueprint(not (y > x)) script_18_6.out [ah@hobbes Documents]\$ python script_18_6.py FalseTrueTrue[ah@hobbes Documents]\$

## Script 18.7

 script_18_7.py # Script 18.7x = 'Hello'y = "world"print(x,y)  # Print both values to screen script_18_7.out [ah@hobbes Documents]\$ python script_18_7.py ('Hello', 'world')[ah@hobbes Documents]\$

## Script 18.8

 script_18_8.py # Script 18.8x = '''This text flows from one lineto the next line inside triple quotes'''print(x) script_18_8.out [ah@hobbes Documents]\$ python script_18_8.py This text flows from one lineto the next line inside triple quotes[ah@hobbes Documents]\$

## Script 18.9

 script_18_9.py # Script 18.9x = 'abcde'print(len(x)) # 5 print(x[0])   # First character 'a' at index 0print(x[-1])  # Last character 'e'print(x[1:3]) # 'bcd' range from index 1, up to, but not including 3 script_18_9.out [ah@hobbes Documents]\$ python script_18_9.py 5aebc[ah@hobbes Documents]\$

## Script 18.10

 script_18_10.py # Script 18.10name = 'Mary'x = 0.12y = 34121.0a = 'Hello {}'.format(name)    # name replaces bracketsprint(a)                       # 'Hello Mary'b = 'X {:.5f} Y {:.2e}'.format(x, y)   # 5 dp float 2 dp sciprint(b)                               # 'X 0.12000 Y 3.41e4' c = '{:3d},{:3d}'.format(12,7)   # pad digits to 3 charactersprint(c)                         # ' 12,  7' script_18_10.out [ah@hobbes Documents]\$ python script_18_10.py Hello MaryX 0.12000 Y 3.41e+04 12,  7[ah@hobbes Documents]\$

## Script 18.11

 script_18_11.py # Script 18.11x = ['a', 'b', 'c', 'd'] print(x[2])   # 'c' - the character at index 2print(x[1:])  # ['b', 'c', 'd'] - from index 1 to the endx[0] = 'z'    # character at index 0 is set to 'z'del x[2]      # Deletes item at index 2print(x)      # ['z', 'b', 'd']y = [[4,7], [9,6]] # A list containing two other listsprint(y[1])        # [9,6] - the sub-list at y index 1print(y[1][0])     # 9 - item 0 from sub-list 1 script_18_11.out [ah@hobbes Documents]\$ python script_18_11.py c['b', 'c', 'd']['z', 'b', 'd'][9, 6]9[ah@hobbes Documents]\$

## Script 18.12

 script_18_12.py # Script 18.12x = list('PQRST')     # Create list from text stringprint(x)              # ['P','Q','R','S','T'] - list of charactersy = list(range(10))   # A list of the range from 0 up to <10 print(y)              # [0,1,2,3,4,5,6,7,8,9] script_18_12.out [ah@hobbes Documents]\$ python script_18_12.py ['P', 'Q', 'R', 'S', 'T'][0, 1, 2, 3, 4, 5, 6, 7, 8, 9][ah@hobbes Documents]\$

## Script 18.13

 script_18_13.py # Script 18.13x = [4, 1, 0] # Three items a, b, c = x   # a is 4, b is 1, c is 0x = ['a', 'b', 'c', 'd']print(len(x))     # 4 - number of items in listprint('c' in x)   # True - string 'c' is in the listprint(2 in x)     # False - number 2 is not in list script_18_13.out [ah@hobbes Documents]\$ python script_18_13.py 4TrueFalse[ah@hobbes Documents]\$

## Script 18.14

 script_18_14.py # Script 18.14x = (9, 3, 1, 0)  # Create a tupley = (2,)          # Tuple with one item (needs trailing comma)z = ()            # Empty tupleprint(x[-1])      # Print the last item - 0#x[0] = 6          # Does not work! Tuples cannot be changed!w = list(x)       # Create equivalent list from a tuplew[0] = 6          # List copy can be changed print('x',x)print('w',w) script_18_14.out [ah@hobbes Documents]\$ python script_18_14.py 0('x', (9, 3, 1, 0))('w', [6, 3, 1, 0])[ah@hobbes Documents]\$

## Script 18.15

 script_18_15.py # Script 18.15x = {1,2,3,4,3,2,1}  # x is {1,2,3,4} - duplicates ignoredy = set([3,4,5,6])   # y created from a listprint(len(y))        # 4print(1 in y)        # False; 1 is not in yprint(x & y)         # {3,4} - intersection, items in bothprint(x | y)         # {1,2,3,4,5,6} - union, items in either script_18_15.out [ah@hobbes Documents]\$ python script_18_15.py 4Falseset([3, 4])set([1, 2, 3, 4, 5, 6])[ah@hobbes Documents]\$

## Script 18.16

 script_18_16.py # Script 18.16d = {"G":329.21, "C":289.18, "A":313.21, "T":314.19}print(d['A'])      # 313.21 -  value associated with 'A'print(len(d))      # 4 - number of key:value pairsprint(d.keys())    # Just keys 'G', 'A', 'T', 'C'print(d.values())  # Just values 329.21, 313.21, 314.19, 289.18 script_18_16.out [ah@hobbes Documents]\$ python script_18_16.py 313.214['A', 'C', 'T', 'G'][313.21, 289.18, 314.19, 329.21][ah@hobbes Documents]\$

## Script 18.17

 script_18_17.py # Script 18.17d = {"G":329.21, "C":289.18, "A":313.21, "T":314.19}d['T'] = 304.19   # Change the value of an existing itemd['U'] = 291.08   # Add a new key:value pairprint(len(d))     # 5 - dict is larger  del d['U']        # Delete a key and its value from the dictionary script_18_17.out [ah@hobbes Documents]\$ python script_18_17.py 5[ah@hobbes Documents]\$

## Script 18.18

 script_18_18.py # Script 18.18x = ['Mon', 'Tue', 'Wed'] # A list of stringsy = ['Fri', 'Sat', 'Sun'] # And anotherx.append('Thu')  # Add a single new item to endprint(x)x.extend(y)       # Extend with items from another collection print(x)x.sort()         # Sort contents alphabeticallyprint(x)x.remove('Sun')  # Remove an itemprint(x)x.index('Sat')   # Positional index of an itemprint(x) script_18_18.out [ah@hobbes Documents]\$ python script_18_18.py ['Mon', 'Tue', 'Wed', 'Thu']['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']['Fri', 'Mon', 'Sat', 'Sun', 'Thu', 'Tue', 'Wed']['Fri', 'Mon', 'Sat', 'Thu', 'Tue', 'Wed']['Fri', 'Mon', 'Sat', 'Thu', 'Tue', 'Wed'][ah@hobbes Documents]\$

## Script 18.19

 script_18_19.py # Script 18.19s = {'G', 'C', 'A', 'T'}   # A set with 4 stringst = {'N', 'R', 'Y'}print(s)s.add('U')       # Add a single item (if not present) print(s)s.update(t)      # Add any new items from another collectionprint(s) script_18_19.out [ah@hobbes Documents]\$ python script_18_19.py set(['A', 'C', 'T', 'G'])set(['A', 'C', 'U', 'T', 'G'])set(['A', 'C', 'G', 'N', 'R', 'U', 'T', 'Y'])[ah@hobbes Documents]\$

## Script 18.20

 script_18_20.py # Script 18.20x = 3if x < -1 or x > 1:    # Run the next indented lines only when true       x *= 2                    print('Value was doubled')  print('Value is:', x)  # Always executed, not in indented block script_18_20.out [ah@hobbes Documents]\$ python script_18_20.py Value was doubled('Value is:', 6)[ah@hobbes Documents]\$

## Script 18.21

 script_18_21.py # Script 18.21x = 3print(x)if x > 0:    print('Positive')elif x < 0:            # Checked if first condition was false    print('Negative')else:                  # If all fails    print('Zero')x = -5print(x)if x > 0:    print('Positive')elif x < 0:            # Checked if first condition was false    print('Negative')else:                  # If all fails    print('Zero') script_18_21.out [ah@hobbes Documents]\$ python script_18_21.py 3Positive-5Negative[ah@hobbes Documents]\$

## Script 18.22

 script_18_22.py # Script 18.22x = 'abc'if x:              # Test innate truth of x    print('true')  # This prints; a non-empty string is trueelse:    print('false') script_18_22.out [ah@hobbes Documents]\$ python script_18_22.py true[ah@hobbes Documents]\$

## Script 18.23

 script_18_23.py # Script 18.23total = 0data = [1,4,9,25,36]for x in data:       # x is first 1, then 4, then 9 etc.    print(x)         # Current value of x in this cycle    total += x       # Add current value of x to total script_18_23.out [ah@hobbes Documents]\$ python script_18_23.py 1492536[ah@hobbes Documents]\$

## Script 18.24

 script_18_24.py # Script 18.24text = 'AGCAGTAGACGAACAT'     # String of charactersfor i, x in enumerate(text):  # Extract index and character value    print(i, x)               # Print index and value for each cycle script_18_24.out [ah@hobbes Documents]\$ python script_18_24.py (0, 'A')(1, 'G')(2, 'C')(3, 'A')(4, 'G')(5, 'T')(6, 'A')(7, 'G')(8, 'A')(9, 'C')(10, 'G')(11, 'A')(12, 'A')(13, 'C')(14, 'A')(15, 'T')[ah@hobbes Documents]\$

## Script 18.25

 script_18_25.py # Script 18.25x = 1while x < 1000:   # Repeat the indented block while this is true    print(x)     x *= 2        # Double the valueprint(x)  # 1024 - final value stopped the loop: not less than 1000 script_18_25.out [ah@hobbes Documents]\$ python script_18_25.py 12481632641282565121024[ah@hobbes Documents]\$

## Script 18.26

 script_18_26.py # Script 18.26t = 0data = [3, -1, 2, -5, None, 9, -2]for x in data:     if x < 0:       continue      # Skip the remainder of 'for' loop   elif x is None:       break         # Quit entirely        t += x * x        # Otherwise do a calculation   print(t)print('Final',t) script_18_26.out [ah@hobbes Documents]\$ python script_18_26.py 91394('Final', 94)[ah@hobbes Documents]\$

## Script 18.27

 script_18_27.py # Script 18.27squares = []for x in range(1, 10):    squares.append(x*x)  # 1, 4, 9, 16, 25, 36...    print('Option 1',squares)squares2 = []squares2 = [x*x for x in range(1,10)]print('Option 2',squares2) script_18_27.out [ah@hobbes Documents]\$ python script_18_27.py ('Option 1', [1])('Option 1', [1, 4])('Option 1', [1, 4, 9])('Option 1', [1, 4, 9, 16])('Option 1', [1, 4, 9, 16, 25])('Option 1', [1, 4, 9, 16, 25, 36])('Option 1', [1, 4, 9, 16, 25, 36, 49])('Option 1', [1, 4, 9, 16, 25, 36, 49, 64])('Option 1', [1, 4, 9, 16, 25, 36, 49, 64, 81])('Option 2', [1, 4, 9, 16, 25, 36, 49, 64, 81])[ah@hobbes Documents]\$

## Script 18.28

 script_18_28.py # Script 18.28x = 1y = 0try:                 # Run the following block and check for failure    w = x / yexcept ZeroDivisionError as err:        # y was zero   print('divided by zero, continuing') # warn, but otherwise ignoreexcept Exception as err:   # Other, unspecified problem    raise(err)             # Trigger the error, do not continue script_18_28.out [ah@hobbes Documents]\$ python script_18_28.py divided by zero, continuing[ah@hobbes Documents]\$

## Script 18.29

 script_18_29.py # Script 18.29def my_calc(x, y=1):     # x is mandatory, y defaults to 1 if not given    z = x * x - y * y    return za = my_calc(7)           # x is 7, y defaults to 1  b = my_calc(x=2, y=2)    # name both argumentsc = my_calc(y=9, x=-1)   # name arguments in a different orderd = my_calc(3, y=-2)     # unnamed arguments (x is 3) come firstprint(a, b, c, d) script_18_29.out [ah@hobbes Documents]\$ python script_18_29.py (48, 0, -80, 5)[ah@hobbes Documents]\$

## Script 18.30

 script_18_30.py # Script 18.30class Person():         # Next indented block is in the class definition  def __init__(self, name, age):  # Values specified when object is made    self.name = name              # Link input values to the object    self.age = age  def get_first_name(self):       # A second, custom function    names = self.name.split()     # self refers to the run-time object    return names[0]               # Give back first wordp1 = Person('Lisa Simpson', 8)    # Make object of Person classp2 = Person('Bart Simpson', 10)   # Make anotherprint(p1.age, p2.age)             # Values linked to objects - 8, 10 print(p1.get_first_name())        # Run a linked function - gives 'Lisa' script_18_30.out [ah@hobbes Documents]\$ python script_18_30.py (8, 10)Lisa[ah@hobbes Documents]\$

## Script 18.31

 script_18_31.py # Script 18.31import math           # Import the inbuilt mathematics moduleprint(math.e)         # An attribute representing eprint(math.exp(2.0))  # Use the exponent (ex) function script_18_31.out [ah@hobbes Documents]\$ python script_18_31.py 2.718281828467.38905609893[ah@hobbes Documents]\$

## Script 18.32

 script_18_32.py # Script 18.32import math as m      # Import as a different nameprint(m.log(2.0))     # Use the logarithm function script_18_32.out [ah@hobbes Documents]\$ python script_18_32.py 0.69314718056[ah@hobbes Documents]\$

## Script 18.33

 script_18_33.py # Script 18.33from math import sqrt, cos   # Import named module componentsx = sqrt(3/2)                # Use the square root functionprint(cos(x))                # Use the cosine function script_18_33.out [ah@hobbes Documents]\$ python script_18_33.py 0.540302305868[ah@hobbes Documents]\$

## Script 18.34

 script_18_34.py # Script 18.34import sysprint(sys.argv) # List of words typed at the command line after "python"print(sys.path) # Directory search path used to locate Python modules sys.exit()      # Quit the Python program script_18_34.out [ah@hobbes Documents]\$ python script_18_34.py ['script_18_34.py']['/home/ah/Documents', '/usr/lib64/python27.zip', '/usr/lib64/python2.7', '/usr/lib64/python2.7/plat-linux2', '/usr/lib64/python2.7/lib-tk', '/usr/lib64/python2.7/lib-old', '/usr/lib64/python2.7/lib-dynload', '/usr/lib64/python2.7/site-packages', '/usr/lib64/python2.7/site-packages/Numeric', '/usr/lib64/python2.7/site-packages/gtk-2.0', '/usr/lib64/python2.7/site-packages/wx-3.0-gtk3', '/usr/lib/python2.7/site-packages'][ah@hobbes Documents]\$

## Script 18.35

 script_18_35.py # Script 18.35import randomrandom.seed(3)                   # Set the random number seedprint(random.randint(1,10))      # A random integer from 1 to 10 inc.print(random.uniform(0.0, 2.0))  # A random float between 0.0 and 2.0data = [1,2,3,4,5]random.shuffle(data)             # Shuffle list order script_18_35.out [ah@hobbes Documents]\$ python script_18_35.py 31.08845845059[ah@hobbes Documents]\$

## Script 18.36

 script_18_36.py # Script 18.36import repattern = re.compile('\d+') # Make a regular expression object                            # (one or more digits)text = 'A 123 B 456'        # A string to look in match_obj = pattern.search(text)  # Match inside string  print(match_obj.group())          # '123' - matching substringprint(match_obj.start)            # 2 - position of matchtext_2 = pattern.sub(text, '**')  # Substitute all matches with '**'print(text_2)                     # 'A ** B **' hits = pattern.findall(text)      # List of matching sub-stringsprint(hits)                       # ['123', '456']) script_18_36.out [ah@hobbes Documents]\$ python script_18_36.py 123**['123', '456'][ah@hobbes Documents]\$

## Script 18.37

 script_18_37.py # Script 18.37import sysprint(sys.path)                        # Current module search pathsys.path.append('./')                  # Directory with my_module.pyimport my_module                       # Import custom module (no .py)print(my_module.CONSTANT)       # Use variable from other Python fileprint(my_module.my_calc(9,-1))  # Use function from other Python file script_18_37.out [ah@hobbes Documents]\$ python script_18_37.py ['/home/ah/Documents', '/usr/lib64/python27.zip', '/usr/lib64/python2.7', '/usr/lib64/python2.7/plat-linux2', '/usr/lib64/python2.7/lib-tk', '/usr/lib64/python2.7/lib-old', '/usr/lib64/python2.7/lib-dynload', '/usr/lib64/python2.7/site-packages', '/usr/lib64/python2.7/site-packages/Numeric', '/usr/lib64/python2.7/site-packages/gtk-2.0', '/usr/lib64/python2.7/site-packages/wx-3.0-gtk3', '/usr/lib/python2.7/site-packages']1.0545718e-3498[ah@hobbes Documents]\$

## Script 18.38

 script_18_38.py # Script 18.38def my_calc(x, y):             # Example function    return x*x - 2*x*y - y*yif __name__ == '__main__':     # Is the code run as a main program?    print(my_calc(2,3))        # Run test code only when it is main    print(my_calc(-1,0))       # Module imports do not run this block script_18_38.out [ah@hobbes Documents]\$ python script_18_38.py -171[ah@hobbes Documents]\$

## Script 18.39

 script_18_39.py # Script 18.39path = 'data_file.txt'            # Location of data in file systemfile_obj = open(path, 'r')        # Create file object for readingdata = file_obj.read()            # Read all the data (as a string)print(data) script_18_39.out [ah@hobbes Documents]\$ python script_18_39.py Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.[ah@hobbes Documents]\$

## Script 18.40

 script_18_40.py # Script 18.1mass = 3.4volume = 1.8density = mass/volume       # Divisiondensity = round(density, 3) # Round to three decimal placesprint(density) script_18_40.out [ah@hobbes Documents]\$ python script_18_40.py ** First read **Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.** Second read **** Third read **Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.[ah@hobbes Documents]\$

## Script 18.41

 script_18_41.py # Script 18.41path = 'data.txt'           # Location of data in file systemfile_obj = open(path, 'r')  # Create file object for readingfor line in file_obj:       # Loop through all lines in file    print(line) script_18_41.out [ah@hobbes Documents]\$ python script_18_41.py chr        mb_size        proteins1        249.250        20122        243.199        12033        198.022        1040[ah@hobbes Documents]\$

## Script 18.42

 script_18_42.py # Script 18.42f_name = 'data.txt'total = 0with open(f_name) as file_obj:    file_obj.readline()         # Read first header line; not used    for line in file_obj:       # Go through each remaining line        data = line.split()     # Split line at whitespace into a list        size = float(data[1])   # Make a number from column [1]        total += size           # Increase total        print(size,total) script_18_42.out [ah@hobbes Documents]\$ python script_18_42.py (249.25, 249.25)(243.199, 492.449)(198.022, 690.471)[ah@hobbes Documents]\$

## Script 18.43

 script_18_43.py # Script 18.43import sys                 # Access the sys modulepy_script = sys.argv[0]    # 'script_18_43.py'   data_file = sys.argv[1]    # 'inputFile.txt'data = open(data_file).read()  # Make file object and read its dataprint(data_file)print(data) script_18_43.out [ah@hobbes Documents]\$ python script_18_43.py inputFile.txt inputFile.txt  SampleID Gender Weight Age Affected1       S1   Male   66.3  20     TRUE2       S2   Male   80.2  21    FALSE3       S3 Female   53.1  29     TRUE4       S4 Female   21.3  18     TRUE5       S5 Female   54.6  NA    FALSE6       S6 Female   78.1  23     TRUE[ah@hobbes Documents]\$

## Script 18.44

 script_18_44.py # Script 18.44from Bio import SeqIO    # Load BioPython module; must be installed file_name = "demoSequences.fasta"     # Location of datafile_obj = open(file_name, "rU")               # Universal read modefor protein in SeqIO.parse(file_obj, 'fasta'):  # Go through each entry  print(protein.id)                            # The ID of seq record  print(protein.seq)                           # The actual sequencefile_obj.close() script_18_44.out [ah@hobbes Documents]\$ python script_18_44.py sp|P52907|CAPZA1_HUMANMADFDDRVSDEEKVRIAAKFITHAPPGEFNEVFNDVRLLLNNDNLLREGAAHAFAQYNMDQFTPVKIEGYEDQVLITEHGDLGNSRFLDPRNKISFKFDHLRKEASDPQPEEADGGLKSWRESCDSALRAYVKDHYSNGFCTVYAKTIDGQQTIIACIESHQFQPKNFWNGRWRSEWKFTITPPTAQVVGVLKIQVHYYEDGNVQLVSHKDVQDSLTVSNEAQTAKEFIKIIENAENEYQTAISENYQTMSDTTFKALRRQLPVTRTKIDWNKILSYKIGKEMQNAsp|Q9Y3C0|CCD53_HUMANMDEDGLPLMGSGIDLTKVPAIQQKRTVAFLNQFVVHTVQFLNRFSTVCEEKLADLSLRIQQIETTLNILDAKLSSIPGLDDVTVEVSPLNVTSVTNGAHPEATSEQPQQNSTQDSGLQESEVSAENILTVAKDPRYARYLKMVQVGVPVMAIRNKMISEGLDPDLLERPDAPVPDGESEKTVEESSDSESSFSD1V1A.pdb_AMLEVVTAGEPLVALVPQEPGHLRGKRLLEVYVGGAEVNVAVALARLGVKVGFVGRVGEDELGAMVEERLRAEGVDLTHFRRAPGFTGLYLREYLPLGQGRVFYYRKGSAGSALAPGAFDPDYLEGVRFLHLSGITPALSPEARAFSLWAMEEAKRRGVRVSLDVNYRQTLWSPEEARGFLERALPGVDLLFLSEEEAELLFGRVEEALRALSAPEVVLKRGAKG-AWAFVDGRRVEGSAFAVEAVDPVGAGDAFAAGYLAGAVWGLPVEERLRLANLLGASVAASRGDHEGAPYREDLEVLLKH.MGKDNYDVVVIGAVGVDTNVYLPGQDIDFDVEANFTTNIDYPGLAGSYTSRGFARLVDKVAVIAYIGEDHNGWYVKNELEGDGIDCTFFLDPVGTKRSINFMYKDGRRKNFYDGKASMEVKPDIDICKSILKKTNLAHFNIVNWTRYLLPVARELGVTISCDIQDITDLDDEYRQDFINYADVLFFSAVNFKDPVPLIKEFIKRKPDRIVICGMGDKGCALGTEQGIKFYKALDMLEPVVDTNGAGDGLAVGFLSSYFIDGFSLEDSILRGQIVARYTCTKKATSSELITGDKLNKYFEEMKESG[ah@hobbes Documents]\$

## Script 18.45

 script_18_45.py # Script 18.45import pandasdata_set = pandas.read_csv('in.csv', sep='\t',                           header=None, names=['A','B','C'])for col in data_set:            # Loop through column names  for value in data_set[col]:   # Column names act as keys    print(col, value) script_18_45.out [ah@hobbes Documents]\$ python script_18_45.py ('A', 'chr')('A', '1')('A', '2')('A', '3')('B', 'mb_size')('B', '249.250')('B', '243.199')('B', '198.022')('C', 'proteins')('C', '2012')('C', '1203')('C', '1040')[ah@hobbes Documents]\$

## Script 18.46

 script_18_46.py # Script 18.46import sqlite3                  # Installed as standardconnection = sqlite3.connect('genome_example.sqlite')cursor = connection.cursor()    # Make a cursormake_table_smt = """ CREATE TABLE Genome (   build_id VARCHAR(12),  genus TEXT NOT NULL,   species TEXT NOT NULL,  strain TEXT,  num_chromsomes INT,  mb_length FLOAT,  pubmed_id VARCHAR(12),  PRIMARY KEY (build_id));"""cursor.execute(make_table_smt)  # Execute SQL string to make tablecursor.close()                  # The cursor is finishedconnection.commit()             # Commit the changes to the database script_18_46.out [ah@hobbes Documents]\$ ls -al genome_example.sqlite -rw-r--r--. 1 ah ah 3072 Jan  5 18:19 genome_example.sqlite[ah@hobbes Documents]\$

## Script 18.47

 script_18_47.py # Script 18.47import sqlite3                  # Installed as standardconnection = sqlite3.connect('genome_example.sqlite')cursor = connection.cursor()    # A new cursor into the database# Create the INSERT statement using "\" to split over two linesadd_genome_smt = 'INSERT INTO Genome (build_id, genus, species, ' \                 'strain, num_chromsomes) VALUES (?, ?, ?, ?, ?)' mouse_data = ['GRCm38', 'Mus', 'musculus', 'C57BL/6', 22] # Data to addcursor.execute(add_genome_smt, mouse_data)                # Add one rowmulti_data = [['GRCh38', 'Homo', 'sapiens', None, 25],    # Data to add              ['GRCh37', 'Homo', 'sapiens', None, 25],              ['BDGP6', 'Drosophila', 'melanogaster', '2057', 7]]cursor.executemany(add_genome_smt, multi_data)            # Add many rowsset_size_smt = 'UPDATE Genome SET mb_length=? WHERE build_id=?'cursor.execute(set_size_smt, [142.573, 'BDGP6'])          # Alter a row cursor.execute(set_size_smt, [3482.010, 'GRCm38'])get_organism_smt = 'SELECT genus, species FROM Genome WHERE build_id=?'genus, species = cursor.execute(get_organism_smt, ['GRCh37']).fetchone()smt = 'SELECT build_id, genus, species, mb_length FROM Genome'for build, g, s, sz in cursor.execute(smt):  print(build, g, s, sz)   cursor.close() script_18_47.out [ah@hobbes Documents]\$ python script_18_47.py (u'GRCm38', u'Mus', u'musculus', 3482.01)(u'GRCh38', u'Homo', u'sapiens', None)(u'GRCh37', u'Homo', u'sapiens', None)(u'BDGP6', u'Drosophila', u'melanogaster', 142.573)[ah@hobbes Documents]\$

## Script 18.48

 script_18_48.py # Script 18.48from numpy import array       # Import the function to make ndarrays x = array([[1,2,3],[4,5,6]])  # Make array from list of listsprint(x.shape)                # (2,3) - rows x columnsprint(x.size)                 # 6 - elements in totalprint(x.ndim)                 # 2 - two axesprint(x.dtype)                # int - whole numbersy = array(x)                  # Copy an arrayy.reshape((3,2))              # [[1,2], [3,4], [5,6]] - Convert to 3 x 2 z = array((3,5,1,2), float)   # Make array forced as floating pointprint('y',y)print('z',z) script_18_48.out [ah@hobbes Documents]\$ python script_18_48.py (2, 3)62int64('y', array([[1, 2, 3],       [4, 5, 6]]))('z', array([ 3.,  5.,  1.,  2.]))[ah@hobbes Documents]\$

## Script 18.49

 script_18_49.py # Script 18.49from numpy import array # Import the function to make ndarrays x = array([5.1, 6.2, 7.3, 8.4, 9.5]) # Flat array of five float numbersx[2]         # 7.3                    - number at index 2x[1:4]       # array([6.2, 7.3, 8.4]) - slice range makes new arrayprint('x',x)y = array([[1,2,3], [4,5,6]]) # Two dimensional arrayy[1,2]       # 6                       - row one, column twoy[0]         # array([1, 2, 3])        - entire row zeroy[0,:]       # array([1, 2, 3])        - row zero, all columns (same)y[:,2]       # array([3, 6])           - all rows, column twoy[-1,:]      # array([4, 5, 6])        - last row, all columnsy[:,1:]      # array([[2, 3],[5, 6]])  - column one onwards y[1,0:2]     # array([4, 5])           - row one, first two columnsy[::-1,:]    # array([[4, 5, 6],[1, 2, 3]]) - reversed rowsy[:,(2,1,0)] # array([[3, 1, 2],[6, 4, 5]]) - new column orderprint('y',y) script_18_49.out [ah@hobbes Documents]\$ python script_18_49.py ('x', array([ 5.1,  6.2,  7.3,  8.4,  9.5]))('y', array([[1, 2, 3],       [4, 5, 6]]))[ah@hobbes Documents]\$

## Script 18.50

 script_18_50.py # Script 18.50import numpy as npa = np.zeros((2,3))          # 2 x 3 array full of 0.0b = np.ones((3,2), int)      # 3 x 2 array full of 1c = np.full((2,2), 7.0)      # 2 x 2 array full of 7.0d = np.identity(3)           # 3 x 3 identity (1 on diagonal 0 elsewhere) e = np.arange(1.0, 3.0, 0.5) # [1.0, 1.5, 2.0, 2.5]                             # From 1.0 up to <3.0 in steps of 0.5print('a',a)print('b',b)print('c',c)print('d',d)print('e',e)x = np.array([1.0, 2.0, 3.0])y = np.array([3.0, 4.0, 5.0])x + 2    # array([3.0, 4.0, 5.0]) - Add 2 to all valuesy / 2    # array([1.5, 2.0, 2.5]) - Divide all values by 2x + y    # array([4.0, 6.0, 8.0]) - i.e. 1+3, 2+4, 3+5 x * y    # array([3.0, 8.0, 15.0])    x - y    # array([-2.0, -2.0, -2.0])x / y    # array([0.33333333, 0.5, 0.6])angles  = np.array([30.0, 60.0, 90.0, 135.0])radians = np.radians(angles)cosines = np.cos(radians)     # array([0.866, 0.50, 0.0, -0.707])print('angles',angles)print('radians',radians)print('cosines',cosines)np.log([10.0, 2.71828, 1.0])  # array([2.302585, 1.0, 0.0])np.exp([2.302585, 1.0, 0.0])  # array([10.0, 2.71828, 1.0]) script_18_50.out [ah@hobbes Documents]\$ python script_18_50.py ('a', array([[ 0.,  0.,  0.],       [ 0.,  0.,  0.]]))('b', array([[1, 1],       [1, 1],       [1, 1]]))('c', array([[ 7.,  7.],       [ 7.,  7.]]))('d', array([[ 1.,  0.,  0.],       [ 0.,  1.,  0.],       [ 0.,  0.,  1.]]))('e', array([ 1. ,  1.5,  2. ,  2.5]))('angles', array([  30.,   60.,   90.,  135.]))('radians', array([ 0.52359878,  1.04719755,  1.57079633,  2.35619449]))('cosines', array([  8.66025404e-01,   5.00000000e-01,   6.12323400e-17,        -7.07106781e-01]))[ah@hobbes Documents]\$

## Script 18.51

 script_18_51.py # Script 18.51from numpy import arange, exp, fft, random, pil = 0.04w = 0.1t = arange(0.0, 100.0, 1.0)       # A range of time valuesx = exp(2j*pi*w*t) * exp(-l*t)    # Wave equation using complex numbersy = fft.fft(x)                    # Fast Fourier transform array of wavea = random.uniform(0.0, 5.0, 100) # Uniform sample 100 values in [0, 5)b = random.normal(0.0, 2.5, 100)  # Sample normal distribution                                  # Mean 0.0, standard deviation 2.5print(x[0],y[0])print(a[0],b[0]) script_18_51.out [ah@hobbes Documents]\$ python script_18_51.py ((1+0j), (0.59324391476340854+1.5043545361088564j))(1.7221136377293216, 3.4608000427150243)[ah@hobbes Documents]\$

## Script 18.52

 script_18_52.py # Script 18.52from matplotlib import pyplotfrom numpy import arange, random  x_vals = arange(1000)y_vals = random.poisson(5, 1000)pyplot.scatter(x_vals, y_vals, s=4, marker='*') # Make scatter plotpyplot.savefig("script_18_52.png", dpi=300)     # Save graph as an imagepyplot.show()                                   # Show on screen script_18_52.out [ah@hobbes Documents]\$ python script_18_52.py [ah@hobbes Documents]\$ ls -al script_18_52.png -rw-rw-r--. 1 ah ah 78409 Jan  5 22:25 script_18_52.png[ah@hobbes Documents]\$ script_18_52.png

## Script 18.53

 script_18_53.py # Script 18.53from scipy import ndimageimg_file = 'my_image.png'pixmap = ndimage.imread(img_file)     # Read image data as arrayheight, width, channels = pixmap.shape           red_channel   = pixmap[:,:,0]         # Color channels are last axisgreen_channel = pixmap[:,:,1]print(red_channel) script_18_53.out [ah@hobbes Documents]\$ python script_18_53.py [[ 0  0  0 ...,  0  0  0] [ 0  0  0 ...,  0  0  0] [ 3  2 14 ...,  8 17 23] ...,  [18 18  7 ..., 18  6  1] [10  9  9 ..., 19  9  7] [11  1  1 ...,  8  1  3]][ah@hobbes Documents]\$

## Script 18.54

 script_18_54.py # Script 18.54def rev_complement(seq):       # Define function name and input argument  rev_map = {'A':'T','T':'A',  # DNA letter mapping dictionary             'G':'C','C':'G'}  rc_seq = ''                  # Rev. comp. string is initially empty  rev_comp = ''  for x in seq:                # Go through each seq letter in turn    y = rev_map[x]             # look up comp. letter in dictionary     rev_comp = y + rev_comp    # Rev. comp. string grows                               # New letter added to _beginning_    return rev_comp              # Rev. comp string sent as output f_seq = 'AGCATAAGAATAGCAGCAGCGCGA'  # A test DNA one-letter sequencer_seq = rev_complement(f_seq)       # Run function on test sequenceprint(f_seq, r_seq)                 # Print original and rev. complement script_18_54.out [ah@hobbes Documents]\$ python script_18_54.py ('AGCATAAGAATAGCAGCAGCGCGA', 'TCGCGCTGCTGCTATTCTTATGCT')[ah@hobbes Documents]\$

## Script 18.55

 script_18_55.py # Script 18.55def calc_seq_ident(seq_a, seq_b):  # Define func. taking two input args  n = min(len(seq_a), len(seq_b))  # Find minimum input seq. length  count = 0.0                      # Starting identity count is zero  for i in range(n):          # Loop through valid position indices    a = seq_a[i]              # Letter at index i for first seq    b = seq_b[i]              # Letter at index i for second seq    if a == b and a != '-':   # Test if letters are same and not a gap      count += 1.0            # .. if they are increase count by one  return 100.0 * count/n      # Calc. and send back identity as % total   seq1 = 'ALIGDPVENTS'seq2 = 'ALIGN-MENTS'x = calc_seq_ident(seq1, seq2) # Run on some test sequencesprint('Seq identity:', x, '%') # x is 72.7 script_18_55.out [ah@hobbes Documents]\$ python script_18_55.py ('Seq identity:', 72.72727272727273, '%')[ah@hobbes Documents]\$

## Script 18.57

 script_18_57.py # Script 18.57def download_db_id(db_id, url, file_name=None): # Save file is optional  #from urllib.request import urlopen    # python 3.0 Import urlopen() function    from urllib2 import urlopen            # python 2.7 Import urlopen() function    response = urlopen(url.format(db_id))  # Makes object to handle link  data = response.read()                 # Read data from URL  data = data.decode('utf-8')            # Interpret data as plain text    if file_name:                       # If save file was specified...    file_obj = open(file_name, 'w')   # Create file object for writing    file_obj.write(data)              # Save all data to file    file_obj.close()  return data                         # Hand back data from functiondef get_uniprot_seq_record(db_id):   import os   from Bio import SeqIO   UNIPROT_FASTA_URL = 'http://www.uniprot.org/uniprot/{}.fasta'   file_name = db_id + '.fasta'  # Add extension to ID to make file name      if not os.path.exists(file_name):     download_db_id(db_id, UNIPROT_FASTA_URL, file_name)      file_obj = open(file_name)    # Open file for reading      for seq_record in SeqIO.parse(file_obj, 'fasta'):     return seq_record           # Give back first record encounteredsr = get_uniprot_seq_record('P18754')print(sr.id, sr.seq) script_18_57.out [ah@hobbes Documents]\$ python script_18_57.py ('sp|P18754|RCC1_HUMAN', Seq('MSPKRIAKRRSPPADAIPKSKKVKVSHRSHSTEPGLVLTLGQGDVGQLGLGENV...EQS', SingleLetterAlphabet()))[ah@hobbes Documents]\$ ls -al P18754.fasta -rw-rw-r--. 1 ah ah 522 Jan  5 22:48 P18754.fasta[ah@hobbes Documents]\$

## Script 18.60

 script_18_60.py # Script 18.60from numpy import random, dot, cov, linalg, concatenate, array from matplotlib import pyplotdef get_principle_components(data, n=2):  data = array(data)               # Convert input to array  mean = data.mean(axis=0)         # Mean vector  centred_data = (data - mean).T   # Centre all vectors and transpose  covar = cov(centred_data)        # Get covariance matrix  evals, evecs = linalg.eig(covar) # Get Eigenvalues and Eigenvectors  indices = evals.argsort()[::-1]  # E. value indices by decreasing size  evecs = evecs[:,indices]    # Sort Eigenvecs according to Eigenvals  p_comp_mat = evecs[:,:n]    # Select required principle components  return p_comp_matsize = (100, 3)                     # 100 points times 3 dimensionsd1 = random.normal(0.0, 0.5, size) d2 = random.normal(0.0, 0.5, size) + array([4.0, 1.0, 2.0])d3 = random.normal(0.0, 0.5, size) + array([2.0, 0.0, -1.0])test_data = concatenate([d1, d2, d3], axis=0)pcomps = get_principle_components (test_data, n=2)   # Run PCApc1, pc2 = pcomps.T   # Extract the two PC vectors from columnsfor x, y, z in (pc1, pc2):  x *= 10             # Scale value so it can be seen better on graph         y *= 10  pyplot.plot((0, x), (0, y))  # Plot PC x, y as line from origin (0,0)x, y, x = test_data.T          # Extract x and y vals from transposepyplot.scatter(x, y, s=20, c='#0080FF', marker='o', alpha=0.5)pyplot.show() script_18_60.out [ah@hobbes Documents]\$ python script_18_60.py [ah@hobbes Documents]\$ script_18_60.png

## Script 18.61

 script_18_61.py # Script 18.61from scipy.stats import norm, poissonfrom numpy import arangefrom matplotlib import pyplotpoisson_rand_var = poisson(2.0)       # Random variable objectx_xals = arange(0, 10, 1)             # Value to plot probabilities fory_vals = poisson_rand_var.pmf(x_xals) # Get probabilities from distrib.pyplot.plot(x_xals, y_vals, color='black') # Make line plot pyplot.show()                              # Show on screen script_18_61.out [ah@hobbes Documents]\$ python script_18_61.py [ah@hobbes Documents]\$ script_18_61.png

## Script 18.62

 script_18_62.py # Script 18.62from numpy import arrayfrom scipy.stats import normdef normal_tail_test(values, mv, std, one_sided=True):    norm_rv = norm(mv, std)     # Normal distrib. random variable object    diffs = abs(values-mv)         # Calc differences from distrib. mean  result = norm_rv.cdf(mv-diffs) # Use cumulative density function    if not one_sided:              # Two-tailed test    result *= 2                  # Distrib. is symmetric: double area  return resultmean    = 1.76std_dev = 0.075values  = array([1.8, 1.9, 2.0])result = normal_tail_test (values, mean, std_dev, one_sided=True)print('Normal one-tail', result)         # [0.297, 0.03097, 0.000687] script_18_62.out [ah@hobbes Documents]\$ python script_18_62.py ('Normal one-tail', array([ 0.29690143,  0.03097408,  0.00068714]))[ah@hobbes Documents]\$

## Script 18.64

 script_18_64.py # Script 18.64def find_discordant_read_pairs(sam_file, max_sep=2e3, verbose=True):    import os  from time import time     # Import inside the function, for a change  from pysam import Samfile, AlignmentFile  start_time = time()       # Record start time (seconds since 1/1/1970)    root_name, ext = os.path.splitext(sam_file) # Split input file ext.  out_file = root_name + "_discordant.bam"    # Out file has new ending    in_sam = Samfile(sam_file, 'rb')          # Open BAM file for reading   out_sam = AlignmentFile(out_file, "wb",   # Make new (empty) BAM file                          template=in_sam)  # base headers on input file      chromo_names = in_sam.references          # Fetch chromosome names  msg_cis = 'Long range: {}:{} - {}:{}'     # Some message templates  msg_trans = 'Interchromosome: {}:{} - {}:{}'    n = 0                    # Number of discordant read pairs found  for align in in_sam:     # Go through each input alignment record    chr_a = align.tid      # Primary chromosome of alignment    pos_a = align.pos      # Primary base pair position    chr_b = align.rnext    # Secondary chromosome    pos_b = align.pnext    # Secondary position    if chr_a == chr_b:          # If both chromosomes the same      delta = abs(pos_b-pos_a)  # Get separation of read ends            if delta > max_sep:         # If too big        out_sam.write(align)      # Write alignment record to output BAM        n += 1                    # Increase count                if verbose:               # Check option to print a message          chr_name = chromo_names[chr_a]          print(msg_cis.format(chr_name, pos_a, chr_name, pos_b))          else:                         # Chromosomes were different      out_sam.write(align)        # Write alignment record to output BAM      n += 1                      # Increase count            if verbose:                 # Check option to print a message        chr_name_a = chromo_names[chr_a]        chr_name_b = chromo_names[chr_b]        print(msg_trans.format(chr_name_a, pos_a, chr_name_b, pos_b))   in_sam.close()  out_sam.close()   delta_time = time()-start_time  # Time taken: current time minus start    # Wite out some useful information  print('Written {:d} records to BAM file: {}'.format(n, out_file))  print('Time taken: {:.1f} s'.format(delta_time)) bam_file = 'GSM409307_UCSD.H3K4me1.bam'  # Test on a _paired_ BAM filefind_discordant_read_pairs(bam_file, verbose=False) script_18_64.out [ah@hobbes Documents]\$ python script_18_64.py [W::bam_hdr_read] EOF marker is absent. The input is probably truncated.Written 6146231 records to BAM file: GSM409307_UCSD.H3K4me1_discordant.bamTime taken: 48.1 s[ah@hobbes Documents]\$ ls -al *.bam-rw-rw-r--. 1 ah ah 356121115 Jan  8 08:04 GSM409307_UCSD.H3K4me1.bam-rw-rw-r--. 1 ah ah 300195982 Jan  8 08:15 GSM409307_UCSD.H3K4me1_discordant.bam[ah@hobbes Documents]\$

## Script 18.65

Not already registered? Create an account now. ×