紧接linux基础(一)的内容
11、安装 miniconda,建立虚拟环境,在虚拟环境中安装生信必备软件
# miniconda官网下载linux版安装包:Miniconda2-4.7.10-Linux-x86_64.sh*
>./Miniconda2-4.7.10-Linux-x86_64.sh* #进入安装包所在目录,运行安装包
>conda config --add channels r
>conda config --add channels conda-forge
>conda config --add channels bioconda
>conda create -n rna python=2 #建立python2环境的名叫rna的虚拟conda环境
>souce activate rna #激活进入rna环境
>conda install -y bwa gatk4 sra-tools fastqc trim-galore star hisat2 bowtie2 subread htseq multiqc samtools salmon #安装软件
>source deactivate #退出环境
12、产生bam文件的命令,示例:
>bowtie2-align-s --wrapper basic-0 -p 20 -x /yanke/yanke21/database/index/bowtie2/hg38 -S WT_rep2_Input.sam -U /tmp/41440.unp
>samtools view tmp.sorted.bam | head -n 2 #运行虚拟环境中的软件要先进入conda环境
SRR1042600.42157053 0 chr1 629895 42 51M * 0 0 ATAACCAATACTACCAATCANTACTCATCATTAATAATCATAATGGCTATA CCCFFFFFHHHHHJJJJJJJ#4AGHJJIIJJIIIIIJJJJIJIIIIJJIJI AS:i:-6 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:11C8A30 YT:Z:UU
SRR1042600.42212881 0 chr1 629895 42 51M * 0 0 ATAACCAATACTACCAATCANTACTCATCATTAATAATCATAATGGCTATA @@<FDFFBFDHHFJEIIGJI#3AFHGEHEIJIIGIIGGIJIIJIGIIGIIJ AS:i:-6 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:11C8A30 YT:Z:UU
SRR1042600.12010763 16 chr1 629895 24 51M * 0 0 ATAACCAATACTTCTAATCAAAACTCATCATTAATAATCATAATGGCTATA ?4B?1*4DD?11*1*?+22+<3F:3@EC:CC4EA,DEDDDDD?D3B:==+; AS:i:-10 XN:i:0 XM:i:4 XO:i:0 XG:i:0 NM:i:4 MD:Z:11C0A1C6T29YT:Z:UU
SRR1042600.29629551 16 chr1 629895 40 51M * 0 0 ATAACCAATACTACCAATCACTACTCATCATTAATAATCATAATGGCTATA HGF?JJHHFDHHGJJIHDFA+E?JIJJIIHGJJJJJJJHHHHHFFFFFCC@ AS:i:-8 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:11C8A30 YT:Z:UU
SRR1042600.41910745 0 chr1 629896 42 51M * 0 0 TAACCAATACTACCAATCAANACTCATCATTAATAATCATAATGGCTATAG CC@FFFFFHHHHGIIHIJJJ#3<CFHCGGIIIJJJJJJJJIGGFHIIJFII AS:i:-6 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:10C9T30 YT:Z:UU
SRR1042600.14329856 16 chr1 629896 8 18M1I32M * 0 0 AAACCAAATCCTCCAATCAAATCCTCATCATTAATAATCATAATGGCTATA #############################@IHHGCE9GHFHHHDDDDD<@@ AS:i:-18 XN:i:0 XM:i:5 XO:i:1 XG:i:1 NM:i:6 MD:Z:0T6T0A2A9A28 YT:Z:UU
13、查看使用的参考基因组 gencode.v29.annotation.gtf.gz 具体有多少条染色体
>zcat gencode.v29.annotation.gtf.gz | grep -v "#" | cut -f 1| sort | uniq -c
250524 chr1
99843 chr10
167540 chr11
159530 chr12
42767 chr13
96124 chr14
100426 chr15
130878 chr16
171137 chr17
48848 chr18
168227 chr19
204338 chr2
61485 chr20
29652 chr21
59764 chr22
172809 chr3
108328 chr4
120515 chr5
124117 chr6
130671 chr7
98583 chr8
97204 chr9
143 chrM
89500 chrX
9064 chrY
14、上面的后缀为BAM 的文件的第二列,只有 0 和 16 两个数字,用 cut/sort/uniq等命令统计它们的个数。
>samtools view tmp.sorted.bam | cut -f 2 | sort | uniq -c
29 0
24 16
15、重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件,再次查看第二列,并且统计
>cd ~/rmDuplicate/samtools/paired
>ll
total 52
drwxrwxr-x 2 yanke21 yanke21 4096 Nov 12 2016 ./
drwxrwxr-x 4 yanke21 yanke21 4096 Nov 12 2016 ../
-rw-rw-r-- 1 yanke21 yanke21 354 Nov 12 2016 readme.txt
-rw-rw-r-- 1 yanke21 yanke21 3111 Nov 12 2016 tmp.header
-rw-rw-r-- 1 yanke21 yanke21 4379 Nov 12 2016 tmp.rmdup.bam
-rw-rw-r-- 1 yanke21 yanke21 1886 Nov 12 2016 tmp.rmdup.vcf.gz
-rw-rw-r-- 1 yanke21 yanke21 11649 Nov 12 2016 tmp.sam
-rw-rw-r-- 1 yanke21 yanke21 4543 Nov 12 2016 tmp.sorted.bam
-rw-rw-r-- 1 yanke21 yanke21 1888 Nov 12 2016 tmp.sorted.vcf.gz
>samtools view tmp.sorted.bam | cut -f 2 | sort | uniq -c
8 147
3 163
1 323
1 353
1 371
1 387
1 433
3 83
2 97
9 99
16、下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件,并且解压,查看里面的文件夹结构, 这个文件有2.3M,注意留心下载时间及下载速度。
>wget -c http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
>unzip sickle-results.zip
>tree sickle-results
sickle-results
├── command.txt
├── single_tmp_fastqc.html
├── single_tmp_fastqc.zip
├── test1_fastqc.html
├── test1_fastqc.zip
├── test2_fastqc.html
├── test2_fastqc.zip
├── trimmed_output_file1_fastqc.html
├── trimmed_output_file1_fastqc.zip
├── trimmed_output_file2_fastqc.html
└── trimmed_output_file2_fastqc.zip
17、解压 sickle-results/single_tmp_fastqc.zip 文件,并且进入解压后的文件夹,找到 fastqc_data.txt 文件,并且搜索该文本文件以 >>开头的有多少行?
>unzip sickle-results/single_tmp_fastqc.zip
>cd single_tmp_fastqc/ && ll
total 296
drwxrwxr-x 4 yanke21 yanke21 4096 Oct 6 2016 ./
drwxr-xr-x 13 yanke21 yanke 4096 Aug 24 11:29 ../
-rw-rw-r-- 1 yanke21 yanke21 12048 Oct 6 2016 fastqc_data.txt
-rw-rw-r-- 1 yanke21 yanke21 5186 Oct 6 2016 fastqc.fo
-rw-rw-r-- 1 yanke21 yanke21 259290 Oct 6 2016 fastqc_report.html
drwxrwxr-x 2 yanke21 yanke21 4096 Oct 6 2016 Icons/
drwxrwxr-x 2 yanke21 yanke21 4096 Oct 6 2016 Images/
-rw-rw-r-- 1 yanke21 yanke21 542 Oct 6 2016 summary.txt
>grep -E '^>>' fastqc_data.txt | wc -l
24
18、下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss文件,去NCBI找到TP53/BRCA1等自己感兴趣的基因对应的 refseq数据库ID,然后找到它们的hg38.tss文件的哪一行。
https://www.ncbi.nlm.nih.gov/gene/7157
>wget -c http://www.biotrainee.com/jmzeng/tmp/hg38.tss
>less hg38.tss
NM_012279 chr5 177020680 177024680 0
NR_131773 chr5 177020680 177024680 0
NR_131774 chr5 177020680 177024680 0
NM_001308213 chr5 177020680 177024680 0
NM_001308219 chr5 177020680 177024680 0
NM_001318082 chr7 99502644 99506644 0
NM_001291980 chr5 177084872 177088872 0
NM_001283039 chr11 3795819 3799819 0
NM_001283040 chr11 3795819 3799819 0
NM_138964 chr2 68641588 68645588 0
NM_139057 chr15 100340005 100344005 1
>grep 'NM_001143990' hg38.tss
NM_001143990 chr17 7684071 7688071 0
>grep 'NM_001291588' hg38.tss
NM_001291588 chr17 43209827 43213827 0
19、解析hg38.tss 文件,统计每条染色体的基因个数。
cut -f 2 hg38.tss | sort | uniq -c
20、解析hg38.tss 文件,统计NM和NR开头的熟练,了解NM和NR开头的含义
grep -E '^NM' hg38.tss | wc -l
51064
grep -E '^NR' hg38.tss | wc -l
15954
NM开头表示转录产物序列,成熟的mRNA序列
NR开头表示非编码的转录子序列,包括RNAs,假基因转子等