Try   HackMD

運用Python解生物資訊問題 (1)

練習一

  1. Calculate the GC content of the following DNA sequence:
    ATGGGTGCGAGAGCGTCAGTATTAAGCGGGGGAGAATTAGATCGATGGGAAAAAATTCGGTTAAGGCCAGGGGGAAAGAAAAAATATAAATTAAAACATATAGTATGGGCAAGCAGGGAGCTAGAACGATTCGCAGTTAATCCTGGCCTGTTAGAAACATCAGAAGGCTGTAGACAAATACTGGGACAGCTACAACCATCCCTTCAGACAGGATCAGAAGAACTTAGATCATTATATAATACAGTAGCAACCCTCTATTGTGTGCATCAAAGGATAGAGATAAAAGACACCAAGGAAGCT
dna_seq = """ATGGGTGCGAGAGCGTCAGTATTAAGCGGGGGAGAATTAGATCGATGGGAAAAAATTCGGTTAAGGCCAGGGGGAAAGAAAAAATATAAATTAAAACATATAGTATGGGCAAGCAGGGAGCTAGAACGATTCGCAGTTAATCCTGGCCTGTTAGAAACATCAGAAGGCTGTAGACAAATACTGGGACAGCTACAACCATCCCTTCAGACAGGATCAGAAGAACTTAGATCATTATATAATACAGTAGCAACCCTCTATTGTGTGCATCAAAGGATAGAGATAAAAGACACCAAGGAAGCT"""
gc = 0
for nt in dna_seq.strip():
    if nt == 'G' or nt == 'C':
        gc += 1
gc_content = gc / len(dna_seq.strip())
print('The GC content of the DNA sequence is %.2f %%' % (gc_content * 100))
  1. Translate the above DNA sequence to protein sequence using the codon table:

codon_table = { 'UUU': 'F', 'UUC': 'F', 'UUA': 'L', 'UUG': 'L', 'UCU': 'S', 'UCC': 'S', 'UCA': 'S', 'UCG': 'S', 'UAU': 'Y', 'UAC': 'Y', 'UAA': '', 'UAG': '', 'UGU': 'C', 'UGC': 'C', 'UGA': '*', 'UGG': 'W', 'CUU': 'L', 'CUC': 'L', 'CUA': 'L', 'CUG': 'L', 'CCU': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P', 'CAU': 'H', 'CAC': 'H', 'CAA': 'Q', 'CAG': 'Q', 'CGU': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R', 'AUU': 'I', 'AUC': 'I', 'AUA': 'I', 'AUG': 'M', 'ACU': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T', 'AAU': 'N', 'AAC': 'N', 'AAA': 'K', 'AAG': 'K', 'AGU': 'S', 'AGC': 'S', 'AGA': 'R', 'AGG': 'R', 'GUU': 'V', 'GUC': 'V', 'GUA': 'V', 'GUG': 'V', 'GCU': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A', 'GAU': 'D', 'GAC': 'D', 'GAA': 'E', 'GAG': 'E', 'GGU': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G', }
Calculate the frequency of each amino-acid, and find the most frequent one.

codon_table = { 'UUU': 'F', 'UUC': 'F', 'UUA': 'L', 'UUG': 'L', 'UCU': 'S', 'UCC': 'S', 'UCA': 'S', 'UCG': 'S', 'UAU': 'Y', 'UAC': 'Y', 'UAA': '', 'UAG': '', 'UGU': 'C', 'UGC': 'C','UGA': '*', 'UGG': 'W', 'CUU': 'L', 'CUC': 'L', 'CUA': 'L', 'CUG': 'L', 'CCU': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P', 'CAU': 'H', 'CAC': 'H', 'CAA': 'Q', 'CAG': 'Q', 'CGU': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R', 'AUU': 'I', 'AUC': 'I', 'AUA': 'I', 'AUG': 'M', 'ACU': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T', 'AAU': 'N', 'AAC': 'N', 'AAA': 'K', 'AAG': 'K', 'AGU': 'S', 'AGC': 'S', 'AGA': 'R', 'AGG': 'R', 'GUU': 'V', 'GUC': 'V', 'GUA': 'V', 'GUG': 'V', 'GCU': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A', 'GAU': 'D', 'GAC': 'D', 'GAA': 'E', 'GAG': 'E', 'GGU': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G', }

transcribe_table = {'A': 'A', 'G': 'G', 'C': 'C', 'T': 'U'}

rna_letters = []
for nt in dna_seq:
    rna_letters.append(transcribe_table[nt])
rna_seq = ''.join(rna_letters)
# print(rna_seq)

protein_seq = []
for i in range(0, len(rna_seq), 3):
    codon = rna_seq[i: i + 3]
    aa = codon_table[codon]
    if aa == '*':
        break
    else:
        protein_seq.append(aa)
protein_seq = ''.join(protein_seq)
print(protein_seq)
  1. Are there other possible start sites in the sequence? How many? For each start site, what would be the sequence of the translated protein? What would be the lengths ofthe DNA and protein sequences? What are the reading frames (0, 1 or 2)?



















start_sites = []
for i in range(len(rna_seq)):
    if rna_seq[i: i + 3] == 'AUG':
        start_sites.append(i)
for start_site in start_sites:
    protein_seq = []
    for i in range(start_site, len(rna_seq), 3):
        codon = rna_seq[i: i + 3]
        if len(codon) == 3:
            aa = codon_table[codon]
        else: 
            break
        if aa == '*':
            break
        else:
            protein_seq += aa
    protein_seq = ''.join(protein_seq)
    report = 'Start site %d: %s, DNA sequence of length %d was translated into %d amino acids; frames = %d' % (start_site, protein_seq, len(rna_seq) - start_site, len(protein_seq), start_site % 3)
    print(report)

練習二

A)
Two strains of bacteria are grown in a lab on a petri dish. Strain A doubles (X2) its population every day,while strain B triples it (X3). The experiment starts with a colony of 10,000 bacteria of type A and 100 of type B on the same petri dish. When the total population (of both strains together) exceeds 108 bacteria,growth is reduced to X1.5 for strain A and X1.8 for strain B. When it exceeds 1010, growth stops. What willbe the number of bacteria during each day of the experiment? How many days will it take until growthstops? What will be the proportion between the two strains at the end of the experiment?




















strain_a = 10000
strain_b = 100
replicate_a = 2
replicate_b = 3
total = strain_a + strain_b
day = 1
while total < 10**10:
    if total > 10**8:
        replicate_a = 1.5
        replicate_b = 1.8
    strain_a *= replicate_a
    strain_b *= replicate_b
    total = strain_a + strain_b
    print('Day%d: strain A: %d bacteria, strain B: %d bacteria [total = %d]' % (day, strain_a, strain_b, total))
    day += 1
prop_a = strain_a / total
prop_b = strain_b / total
print('At Day%d, the total bacteria reached %d with %.2f %% of strain A and %.2f %% of strain B' % (day - 1, total, prop_a * 100, prop_b * 100))

In a second experiment, a bacteriophage is introduced into the dish. When the total population exceeds

106 bacteria, it starts spreading and kills 40% of the cells every day (of both strains). Every day after the bacteriophage becomes active, 10% of type-A bacteria will develop immunity to the bacteriophage and retain normal growth rate. Strain B has a lower mutation rate, so only 1% of its cells will develop immunity every day. Run the simulation again and determine the fraction of the 4 populations at the end of the experiment (type A/B with/without immunity).







































immunity_development_rate = {'A': 0.1, 'B': 0.01}
populations = {('A', False): 1e4, ('A', True): 0, ('B', False): 1e2, ('B', True): 0}
day = 0
is_phage_active = False
while sum(populations.values()) < 1e10:
    if sum(populations.values()) > 1e6:
        is_phage_active = True
    if sum(populations.values()) < 1e8:
        growth_rate = {'A': 2, 'B': 3}
    else: 
        growth_rate = {'A': 1.5, 'B': 1.8}
    for strain_key in populations.keys():
        strain, is_immune = strain_key
        populations[strain_key] *= growth_rate[strain]
        if is_phage_active and not is_immune:
            populations[strain_key] *= 0.6
    if is_phage_active:
        for strain in ['A', 'B']:
            immunity_rate = immunity_development_rate[strain]
            num_developed_immunity = immunity_rate * populations[(strain, False)]
            populations[(strain, True)] += num_developed_immunity
            populations[(strain, False)] -= num_developed_immunity
    day += 1
    if is_phage_active:
        print('After %d days [active phage]:' % day)
    else:
        print('After %d days:' % day)
    for strain in ['A', 'B']:
        for is_immune in [False, True]:
            if is_immune:
                immune_label = 'Immune'
            else:
                immune_label = 'Unimmune'
            population_size = populations[(strain, is_immune)]
            prop = population_size / sum(populations.values()) * 100
            print('\t' + '%s %s: %d bacteria(%.2f%%)' % (immune_label, strain, population_size, prop))
    print('\t' + 'total bacteria: %d' % sum(populations.values()))

練習三

Find the reverse complement of the following DNA sequence and translate the output DNA sequence into protein sequence: CTACTGTAATTCAACACAACTGTTTAATAGTACTTGGTTTAATAGTACTTGGAGTACTGAAGGGTCAAATAACACTGAAGGAAGTGACACAATCACCCTCCCATGCAGAATAAAACAAATTATAAACATGTGGCAGAAAGTAGGAAAAGCAATGTATGCCCCTCCCATCAGTGGACAAATTAGATGTTCATCAAATATTACAGGGCTGCTATTAACAAGAGATGGTGGTAATAGCAACAATGAGTCCGAGATCTTCAGACCTGGAGGAGGAGATATGAGGGACAATTGGAGAAGTGAATTATATAAATATAAAGTAGTAAAAATTGAACCATTAGGAGTAGCACCCACCAAGGCAAAGAGAAGAGTGGTGCAGAGAGAAAAAAGAGCAGTGGGAATAGGAGCTTTGTTCCTTGGGTTCTTGGGAGCAGCAGGAAGCACTATGGGCGCAGCCTCAATGACGCTGACGGTACAGGCCAGACAATTATTGTCTGGTATAGTGCAGCAGCAGAACAATTTGCTGAGGGCTATTGAGGCGCAACAGCATCTGTTGCAACTCACAGTCTGGGGCAT























dna_seq = """CTACTGTAATTCAACACAACTGTTTAATAGTACTTGGTTTAATAGTACTTGGAGTACTGAAGGGTCAAATAACA CTGAAGGAAGTGACACAATCACCCTCCCATGCAGAATAAAACAAATTATAAACATGTGGCAGAAAGTAGGAAAA GCAATGTATGCCCCTCCCATCAGTGGACAAATTAGATGTTCATCAAATATTACAGGGCTGCTATTAACAAGAGA TGGTGGTAATAGCAACAATGAGTCCGAGATCTTCAGACCTGGAGGAGGAGATATGAGGGACAATTGGAGAAGTG AATTATATAAATATAAAGTAGTAAAAATTGAACCATTAGGAGTAGCACCCACCAAGGCAAAGAGAAGAGTGGTG CAGAGAGAAAAAAGAGCAGTGGGAATAGGAGCTTTGTTCCTTGGGTTCTTGGGAGCAGCAGGAAGCACTATGGG CGCAGCCTCAATGACGCTGACGGTACAGGCCAGACAATTATTGTCTGGTATAGTGCAGCAGCAGAACAATTTGC TGAGGGCTATTGAGGCGCAACAGCATCTGTTGCAACTCACAGTCTGGGGCAT"""
dna_seq = dna_seq.replace(" ", "")
reverse_dna_table = {'A': 'T', 'G': 'C', 'T': 'A', 'C': 'G'}
rev_dna_seq = []
for nt in dna_seq[::-1]:
    rev_dna_seq.append(reverse_dna_table[nt])
rev_dna_seq = ''.join(rev_dna_seq)

rna_seq = rev_dna_seq.replace('T', 'U')
protein_seq = []
for i in range(0, len(rna_seq), 3):
    codon = rna_seq[i: i + 3]
    if len(codon) == 3:
        aa = codon_table[codon]
    else:
        break
    if aa == '*':
        break
    else:
        protein_seq.append(aa)
protein_seq = ''.join(protein_seq)
print(f'Reverse DNA sequence: \n{rev_dna_seq} \nProtein sequence: \n{protein_seq}')

What is the frequency of each codon in the obtained sequence?














codon_counts = {}
for i in range(0, len(rna_seq), 3):
    codon = rna_seq[i: i+3]
    if codon in codon_counts:
        codon_counts[codon] += 1
    else:
        codon_counts[codon] = 1

codon_frequencies = {}
total_count = sum(codon_counts.values())
for codon, count in codon_counts.items():
    codon_frequencies[codon] = count / total_count
    print('%s: %.2f' % (codon, codon_frequencies[codon]))

Does the distribution of codons seem uniform? Are there codons that appear dramatically more (or less) than expected (given their amino-acid frequency in that sequence)?



















aa_to_codon = {}
for codon, aa in codon_table.items():
    if aa == '*':
        continue
    if aa in aa_to_codon:
        aa_to_codon[aa].add(codon)
    else:
        aa_to_codon[aa] = {codon}

for aa, aa_codons in aa_to_codon.items():
    aa_counts = protein_seq.count(aa)
    if aa_counts == 0:
        continue
    print('Amino acid %s (%d codons, %d occurrence in the protein sequence' % (aa, len(aa_codons), aa_counts))
    for codon in aa_codons:
        codon_count = codon_counts.get(codon, 0)
        codon_freq = codon_count / aa_counts
        print('\t' + 'Codon %s: %d occurrence (%.1f%%)' % (codon, codon_count, codon_freq * 100))
Last changed by 
 
Atomic Notes
這裡是Eric Juo(小卓)無聊發起的小小空間，目的在於共享生活中學習到的知識及經驗，這裡是一個新手友善空間，不要求筆記寫的多完整，只希望學習到的東西可以在某個地方保存下來，並可以幫助其他人一起學習
0
212

Read more

RNA-Seq overivew slide

lectures and hands on tutorial

Apr 24, 2022
RNA-Seq Analysis

RNA-Seq stands for RNA sequencing. It is a technology for studying RNA species. In a cell, the majority of RNA sepcies is ribosomal RNA (rRNA) which accounts for 90% of total RNA. Conversly, messenger RNA (mRNA) only aocunts for 1-2% of total RNA. The central dogma of biology tells us that DNA transcribes into mRNA and then mRNA translates into protein to execute biological functions. mRNA is the only RNA species we know encodes protein seuqences. Studying the level of mRNAs enables us to understand the gene expression level and infer the protein expression level. Therefore, most of time we only interest in the level of mRNA, and wants to deplet rRNA from the total RNA. Preprocessing RNA samples before sequencing includes total RNA extraction, mRNA enrichment, fragmentation, first-strand cDNA synthesis, RNA strand digestion, second-strand cDNA synthesis, 3&apos;end repair, adenylation and adatper ligation. In total RNA extraction step, extracted RNA should not have signs of degradation. In mRNA enrichment, poly-T beads or poly-T columns are used to separate mRNA from other RNA sepcies, since only mRNA contains poly-A tails. mRNAs are fragmented by chemicals or ultrasonic into certain size range that fits for the sequencing capacity of sequencer. mRNA are primed with random hexamer and reverse-transcribed into first-strand cDNA, followed by RNA strand digestion and second-strand cDNA synthesis. For making strand-specific RNA-seq, TTPs are replaced with dUTPs, which serves as marks for the second-strand, in the second-strand cDNA synthesis step. Since the nature of polymerase synthesis, the double strand cDNA product lose several bases at the 3&apos; ends on each strand. The 3&apos;ends are repaired by dna end repair enzyme and then added a adenine. Y-shape adapters with 3&apos;-T overhang are ligated to the 3&apos;-A overhang of the double strand cDNA. To reduce the complexity in short read assembly, dUTP marked cDNA stands are digested using uracil-N-glycolas (UNG). The single strand cDNAs are then subjected to sequencing. In sequencing step, primer targeting adapter sequence aneal to one end of cDNA. Fluoresence-labelled ATPs, CTPs, GTPs or TTPs is added one at a time and photographed. Layers of photographes are analyzied by computer to interpret which nucleotide is added to the growing strand in each reaction. For paired end sequencing, the other primer targeting adatper on the other side is added and squenced again. The resultant are 2 separate files, one from forward strand, the other from reverse strand. Raw reads obtained from sequencer should undergo quality control before downtream analysis. Raw reads data are normally reported in fastq format. Fastq format stores 4 lines of information for each read. The first line is read name, the second line is sequence, the thrid line is a "+" sign, and the forth line is ASCII characters encypting qualtiy scores. For illumina sequencer, the quality socre is reported in Phred33 scheme. To obtain numberical quality score, ASCII characters need to be converted to numbers they represented in computer and minus 33. Quality score is used to repersent how confidence we can say this base is called correctly. A 40 quality score means 1 in 10000 chances this base is called incorrectly. Thus, higher quality score indicating higher chance this base is called correctely. Usually, we would want to trim off bases with quality score lower than 20. Low quality score base calling usually happens at the 3&apos; end of reads since the polymerase is getting weak attaching to the template and prone to add wrong nucleic acids. FastQC is a commonly used tool to give an overview of read quality. It reports per base quality, per sequence quality, per base nucleotide content, per sequence GC content, number of duplicated sequences, overrepresented sequences and adapter content metrices. Per base quality metric shows the distribution of quality score in box plot for each base location. The upper boundary of the yellow box is the third quantile of qualtiy scores, and the lower boundary of the box is the first quantile. The red line at the middle of box is the median value of quality scores. Per sequence quality metric averages base quality for each read and shows the distribution of per read quality. For good quality reads, the distribution of read quality should skews toward high quality side. Per base nucleotide content metric shows the frequencies of nucleic acid appear at each base location. In principle, RNA has been randomly fragmented, so the frequency of each nucleic acid should distribute evenly across all read length. However, RNA-Seq data is an exception, a research has found that priming with random hexamer in cDNA synthesis step seems to have certain selection for fragmented RNA sequences, which then results in not that random nucleic acid content at the 5&apos;end of reads. But these bases are still real bases in mRNAs, not artifacts, it&apos;s ok to put them into short read assembly. Per sequence GC content metric cacluate GC content for each read and report the distribution of per read GC content. GC content is like a fingerprint of a species, it has a constant value for a given species. This also applys to GC content of reads. The distribution of per read GC content should peaks at the value equal to that species&apos; GC content. If peak shifts away from the expected GC content, it can stem from foreign species contaminants in read dataset. Number of duplicated sequences metric shows how many reads are exact maches to each other. For RNA-seq data set, highly expressed transcripts usually repeatly sequenced resulting in many duplicated reads. Overrepresented sequences metric reports reads that appear repeatively and account for 0.1% of total reads. Adapters are the most likely sequences captured by this metric. Adapter content metric align reads with public adapters and report alignment hits.

Apr 15, 2022
三角函數

Fundamental Trigonometric Identities Power Reducing Formulas Fucdamental Trigonometric Identities Reciprocal identities (&#x5012;&#x6578;) $sec\theta=\frac{1}{cos\theta}$ $csc\theta=\frac{1}{sin\theta}$ $cot\theta=\frac{1}{tan\theta}$ $cos\theta=\frac{1}{sec\theta}$

Oct 11, 2020
FASTAQ format

FastaQ&#x683C;&#x5F0F;&#x662F;&#x7528;&#x4F86;&#x5132;&#x5B58;&#x6838;&#x9178;&#x5E8F;&#x5217;&#x53CA;&#x5176;&#x5C0D;&#x61C9;&#x9E7C;&#x57FA;&#x54C1;&#x8CEA;&#x7684;&#x8CC7;&#x6599;&#x5132;&#x5B58;&#x683C;&#x5F0F;&#x3002; &#x5E8F;&#x5217;&#x54C1;&#x8CEA; &#x5E8F;&#x5217;&#x54C1;&#x8CEA;&#x662F;&#x7528;&#x4F86;&#x8A55;&#x4F30;&#x6B21;&#x4E16;&#x4EE3;&#x5B9A;&#x5E8F;(NGS)&#x4E4B;&#x5B9A;&#x5E8F;&#x7D50;&#x679C;&#x53EF;&#x4FE1;&#x5EA6;&#x7684;&#x6578;&#x503C;&#xFF0C;&#x5176;&#x516C;&#x5F0F;&#x5982;&#x4E0B;: $$Q = -10logP$$ Where Q:&#x4EE3;&#x8868;&#x5E8F;&#x5217;(&#x9E7C;&#x57FA;)&#x54C1;&#x8CEA; P:&#x4EE3;&#x8868;&#x51FA;&#x73FE;&#x932F;&#x8AA4;&#x7684;&#x767E;&#x5206;&#x6BD4;

Aug 23, 2020
Read more from Atomic Notes
Published on HackMD