Posts

Showing posts with the label NGS

What are the best ways to start learning bioinformatics for a wet lab biologist?

Image
I came across this question several times. In my real life, I also encountered this question many times. Many time, my friend asked me: “Hey, can I learn Bioinformatics with your? Can you give me some materials? ” At the beginning, I would say: “Sure. Let me know if you need any help.” Then mostly nothing happened after I sent them the web links or e-books. Now if someone ask me the same question, I usually will say: “Are you sure? If you take it seriously, I’ll teach you.” My point here is take it seriously. What is Bioinformatics? Bioinformatics is the science of collection and analyzing complex biological data such as genetics codes. This is the definition given by the first result in Google search. Maybe it’s too abstract to understand. Based on my understanding, I want to talk about Bioinformatics in the following aspects: Software development Maybe the most popular tool is blast . When I was in college, I selected a course called “生物信息学”. The english name is “Bioinformati...

Software Tips: Deconseq

Problem solved: In the README file, they give a database link:  ftp://edwards.sdsu.edu:7009/ deconseq/db/ It seemed that this link doesn't work. It appears the SDSU firewall is blocking the FTP site.  I moved the databases to an HTTP page for the time being. https://edwards.sdsu.edu/~ katelyn/DECONSEQ/

Software Tips: Deconseq

Problem solved: In the README file, they give a database link:  ftp://edwards.sdsu.edu:7009/ deconseq/db/ It seemed that this link doesn't work. It appears the SDSU firewall is blocking the FTP site.  I moved the databases to an HTTP page for the time being. https://edwards.sdsu.edu/~ katelyn/DECONSEQ/

Shell tips: Can a Bash script tell what directory it's stored in?

Question: How do I get the path of the directory in which a  Bash  script is located FROM that Bash script? Answer: DIR = "$( cd " $ ( dirname "${BASH_SOURCE[0]}" ) " && pwd )" is a useful one-liner which will give you the full directory name of the script no matter where it is being called from. This will work as long as the last component of the path used to find the script is not a symlink (directory links are OK). If you want to also resolve any links to the script itself, you need a multi-line solution: SOURCE = "${BASH_SOURCE[0]}" while [ - h "$SOURCE" ]; do # resolve $SOURCE until the file is no longer a symlink DIR = "$( cd -P " $ ( dirname "$SOURCE" ) " && pwd )" SOURCE = "$(readlink " $SOURCE ")" [[ $SOURCE != /* ]] && SOURCE = "$DIR/$SOURCE" # if $SOURCE was a relative symlink, we need to resolve it relative ...

NGS software tips: SICER for ChIP-seq

Tip: Do NOT need to move SICER anymore. For SICER, if you are going to the original code, you have to copy the source code to your working directory. This is because the PATHTO variable varies when you’re doing different work. To avoid this, you can modify the code to: PATHTO="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" #SICER=$PATHTO/SICER SICER=$PATHTO PYTHONPATH=$SICER/lib export PYTHONPATH Tip: About duplicate reads One step of SICER is to remove duplicate read using remove_redundant_reads.py. “-t 1” is the default cutoff when doing duplicate reads removal. See line 60 of remove_redundant_reads.py . Solved: Error Report.   File "/depot/bioinfo/Projects/Kal/bean_GHvsTC/software/SICER_V1.1/SICER/src/remove_redundant_reads.py", line 84     print chrom, "\tPlus reads:",p_total, "\tRetained plus reads:", p_retained,     ";\tMinus reads:", m_total, "\tRetained...

NGS tips: FPKM for short genes in cufflinks

Don't trust Cufflinks FPKM for short genes Here is what Cole (author of Cufflinks) commented to the observation of " very high RPKM values from Cufflink ": This issue has been discussed elsewhere on this board. As Nicholas points out, RNA-Seq really isn't reliable for very short transcripts. The reason is that all the fragments that map to these transcripts come from the "tail" of the distribution of library fragment lengths. That is, fragments that map to microRNAs are much, much shorter than most fragments in the library - by design in the RNA-Seq protocol, which size selects away very short inserts. Thus,  Cufflinks infers that even though relatively few fragments actually mapped to the microRNAs, there were probably TONS of individual microRNA molecules in the transcriptome before all of the various size selection parts of the protocol kicked in. Cufflinks accordingly increases the FPKM of these short transcripts to compensate for the bias agai...

NGS tips: Picard MarkDuplicates - whole bam file marked as duplicates?

Image
If you used Picard MarkDuplicates to mark duplicate reads in your bam files, you may notice that each line has a tag "PG:Z:MarkDuplicates".  You would be confused. Each read has this mark. Does this mean that every read is duplicate?  The short answer is NO.  In SAM/BAM file, whether a read is duplicate is determined by FLAG (the 2nd column).  You can decode SAM flags under the link: https://broadinstitute.github.io/picard/explain-flags.html For example, I have a read which has a tag "PG:Z:MarkDuplicates". You can type 83 in the box and press "Expain" button. Then you can see the summary shown in the left.  Each read has 12 properties. The 11th is to determine whether this is duplicated or not. As you can see, for 83 the 11th property is not marked meaning that this one is not a duplicate one.  You can also use REMOVE_DUPLICATES = true to remove duplicate reads in the bam files.  ...