常見生信操作

2021-09-12 09:20:10 字數 2127 閱讀 4798

安裝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...