安裝samtools :conda install samtools
# srand: 隨機數發生器。設定固定的種子, 保證每次出來的結果一致
# rand: 返回[0,1)之間的隨機數, 包含0不包含1
1.產生隨機的基因組檔案
echo 1 | awk -v seed=1 -v label=chr -v chrnum=4 -v expected_len=60000 -f generaterandom.awk >genome.fa
2.獲得突變檔案
#check iupac here:
# -n: 獲得40k read pairs
# mut.txt: 突變位點或區域
3.建立基因組索引
bwa index genome.fa
# samtools fadix快速獲取某區域序列
samtools faidx genome.fa
4.序列比對回基因組
bwa mem -t 3 genome.fa ehbio_1.fq ehbio_2.fq | gzip >map.sam.gz
5.篩選比對上的高質量reads
samtools view -f4 -q1 -b map.sam.gz -o map.bam
# 下面2個排序用法都可以, 看使用的samtools版本
samtools sort -@ 2 map.bam map.sortp
samtools sort -@ 2 -o map.sortp.bam map.bam
samtools index map.sortp.bam
6.統計比對reads 數
samtools view -c map.sortp.bam
7.統計未比對上的reads 數
samtools view -c -f 4 map.sam.gz
8.統計比對到正鏈的reads 數
samtools view -c -f 16 map.sam.gz
9.獲取properly-paired 的reads 數
samtools view -f2 -f 256 -c map.sortp.bam
10.檢視每個位置鹼基比對或錯配情況
# -q 0: 測試資料使用, 預設為-q 13, 表示過濾掉低質量測序鹼基
samtools mpileup -f genome.fa -q 0 map.sortp.bam | less
#測序鹼基列解釋:
1. 點(.) 代表匹配正鏈鹼基
2. 逗號(,) 代表匹配負鏈鹼基
3. 大寫字母(acgtn) 表示正鏈錯配
4. 小寫字母(acgtn) 表示負鏈錯配
5. 模式\+[0-9]+[acgtnacgtn]+ 表示在當前參考位置和下乙個參考位置之間有插入,插入鹼基數是+ 後
面的證書,插入鹼基是數字後面的字母串。下面展示的是2 bp 的插入
seq2 156 a 11 .$……+2ag.+2ag.+2aggg <975;:<<<<<
6. 模式-[0-9]+[acgtnacgtn]+' 參考基因組存在鹼基缺失。下面展示的是4 bp『缺失:
seq3 200 a 20 „„,..,.-4cacc.-4cacc….,.„.ˆ~. ==<<<<<<<<<<<::>
7. 符號ˆ 表示測序序列起始位置落於此(a symbol ˆ' marks the start of a read segment which
is a contiguous subsequence on the read separated byn/s/h』 cigar operations). 後面跟
隨的符號的ascii 值減去33 表示該位置鹼基的質量。符號『$』 表示測序序列片段的終止。主要用於從
pileup 檔案中獲得原始測序reads。
安裝bedtools:conda install bedtools
bedtools是處理基因組資訊分析的強大工具集合
重複序列區域的獲取也可以用:
待續。。。。。
生信基礎(二) 生信學習資料
原創 hxj7 上次談到生信人員需要熟練掌握一些程式語言,還講了perl和python的選擇問題。那麼,如果已經選定了一門程式語言,到底該如何學習它呢?今天的我們可以通過mooc跟著名師學習或者上知乎提問,幸運的話還能得到大牛指點。不過,在我剛接觸程式設計的時候,mooc和知乎都還未興起,所以我都是...
生信轉崗心得
應健明的邀請,我也寫一篇關於轉行做生物資訊的心得,本來以為很輕鬆就可以寫出來的,但是發現並不那麼好寫。如果放在去年剛轉崗之時,我想應該更順手,那時感觸良多且深刻。不過我雖是個健忘的人,但是還清醒地記得去年7月份通過公司的轉崗答辯之時,心情無比的愉快與美麗,覺得終於從一名生信愛好者成為一名正規軍。之前...
易生信Linux培訓
生信寶典 和 巨集基因組 聯合舉辦的生物資訊系列9天實戰班,歷時5個月,已經順利完成了4期,包括轉錄組分析 高通量測序分析的基礎課 r語言 cytoscape繪圖 資料視覺化是解釋資料不可缺少的一部分 python程式設計 每個人都應該會一門程式語言 微生物組擴增子測序分析 當前最火的領域之一 li...